diff --git a/Project.toml b/Project.toml index 405463ae..dd1b1b5d 100644 --- a/Project.toml +++ b/Project.toml @@ -25,6 +25,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" NestedUnitRanges = "032820ab-dc03-4b49-91f4-7d58d4da98b3" OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Requires = "ae029012-a4dd-5104-9daa-d747884805df" SauterSchwab3D = "0a13313b-1c00-422e-8263-562364ed9544" @@ -63,8 +64,9 @@ LiftedMaps = "0.5.1" LinearMaps = "3.7 - 3.9, 3.11.2" NestedUnitRanges = "0.2.3" OhMyThreads = "0.8.3" -ProgressMeter = "1.11.0" PlotlyJS = "0.18.18" +Polynomials = "4.1.1" +ProgressMeter = "1.11.0" Requires = "1" SauterSchwab3D = "0.2" SauterSchwabQuadrature = "2.4.0" @@ -74,5 +76,3 @@ Suppressor = "0.2.8" TestItems = "0.1.1, 1" WiltonInts84 = "0.2.8" julia = "1.10" - - diff --git a/docs/make.jl b/docs/make.jl index 093b4867..e4ba3333 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -72,6 +72,9 @@ makedocs(; "Low Frequency Stable PMCHWT"=>"projectors/pmchwt_lf.md", "Dense Grid Stable PMCHWT"=>"projectors/pmchwt_theta.md", ], + "OSRC operators"=>Any[ + "OSRC operators"=>"operators/OSRC_efie.md", + ], ], "Basis Functions" => Any[ "Overview"=>"bases/overview.md", diff --git a/docs/runliterate.jl b/docs/runliterate.jl index 83a1bd9b..fe5cf950 100644 --- a/docs/runliterate.jl +++ b/docs/runliterate.jl @@ -3,4 +3,5 @@ using Literate Literate.markdown("./examples/composedoperator.jl","./docs/src/composedoperator"; execute=true) Literate.markdown("./examples/pmchwt_theta.jl","./docs/src/projectors"; execute=true) Literate.markdown("./examples/pmchwt_lf.jl","./docs/src/projectors"; execute=true) +Literate.markdown("./examples/OSRC_efie.jl","./docs/src/operators"; execute=true) diff --git a/docs/src/operators/OSRC_efie.md b/docs/src/operators/OSRC_efie.md new file mode 100644 index 00000000..cbff12c4 --- /dev/null +++ b/docs/src/operators/OSRC_efie.md @@ -0,0 +1,171 @@ +```@meta +EditURL = "../../../examples/OSRC_efie.jl" +``` + +In this example, the preconditioning ability of the On-Surface Radiation Condition (OSRC) Magnetic-to-Electric (MtE) map is shown on the Electric Field Integral Equation (EFIE), on the sphere. +# OSRC operators +The MtE map ``\text{V}^{+}`` is a map of the trace of the magnetic field ``n \times h`` to the trace of the electric field ``n \times e`` on a scattering surface. +The approximation of the local MtE surface operator as an On-Surface Radiation Condition (OSRC) operator on an arbitrary surface ``\Gamma`` is given by: +```math +\begin{equation} + \text{V}^{+} \left( n \times h \right) + n \times e = 0, \quad \text{on} \ \Gamma, +\end{equation} +``` +In the implementation a highly accurate pseudo-differential operator ``\text{V}_{\epsilon}^{+}`` is used to represent the Magnetic-to-Electric map. +```math +\begin{aligned} + \text{V}_{\epsilon}^{+} &:= - \pmb{\Theta}^{-1} \Lambda_{2, \epsilon}^{-1} \Lambda_{1, \epsilon}, \\ + \pmb{\Theta} \left( \cdot \right) &:= \left( n \times \cdot \right), \\ + \Lambda_{2, \epsilon} &:= I - \textbf{curl}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{curl}_{\Gamma}, \\ + \Lambda_{1, \epsilon} &:= \left(I + \mathcal{J} \right)^{1/2}, \\ + \mathcal{J} &:=\frac{\Delta_{\Gamma}}{\kappa_{\epsilon}^2} = \textbf{Grad}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{Div}_{\Gamma} - \textbf{curl}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{curl}_{\Gamma}. +\end{aligned} +``` +More details on this OSRC operator and its implementation are given in +[An OSRC Preconditioner for the EFIE (Betcke et al. (2021))](https://arxiv.org/abs/2111.10761) +The Electric-to-Magnetic (EtM) OSRC operator is also implemented, as given in +[Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell’s equations (C. Geuzaine et al. (2014))](https://www.researchgate.net/publication/261636307_Approximate_local_magnetic-to-electric_surface_operators_for_time-harmonic_Maxwell's_equations) + +# OSRC EFIE preconditioning example +We begin with loading needed packages + +````julia +using BEAST, CompScienceMeshes + +using Makeitso +using LinearAlgebra + +using Plots +using PlotlyDocumenter +```` + +Now, assemble and solve the unpreconditioned EFIE + +````julia +BLAS.set_num_threads(8) + +@target geo (;h) -> begin + (; Γ = CompScienceMeshes.meshsphere(radius=1.0, h=h)) +end + +@target formulation_EFIE (geo,; κ) -> begin + X = raviartthomas(geo.Γ) + T = Maxwell3D.singlelayer(wavenumber=κ) + + E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ); + e = (n × E) × n; + + bx = assemble(e, X) + A = assemble(T,X,X); + return (;bilforms=(;A), linforms=(;bx)) +end + +@target solution_EFIE (formulation_EFIE,; residual) -> begin + (;bilforms, linforms) = formulation_EFIE + (;A) = bilforms; (;bx) = linforms; + iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500) + u, ch = BEAST.solve(iT, bx) + return (;iters=ch.iters, u=u) +end +```` + +```` +0x4d338de5db5a428f +```` + +Next, assemble the OSRC MtE operator on the primal grid and use it to precondition the EFIE + +````julia +@target OSRC_MtE_op (geo,;κ, Np, curvature) -> begin + MtE_OSRC_operator = BEAST.MtE_OSRC_op(κ, Np, pi/2, curvature) + Nd = BEAST.nedelec(geo.Γ) + MtE_map = assemble(MtE_OSRC_operator, Nd, Nd) + return (;MtE=MtE_map) +end + +@target solution_OSRC_EFIE (formulation_EFIE, OSRC_MtE_op; residual) -> begin + (;bilforms, linforms) = formulation_EFIE + (;A) = bilforms; (;bx) = linforms; + P_OSRC = OSRC_MtE_op.MtE + iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500, left_preconditioner=P_OSRC) + u, ch = BEAST.solve(iT, bx) + return (;iters=ch.iters, u=u) +end +```` + +```` +0x48b1ab5f44e88685 +```` + +Finally, assemble and use the Calderón preconditioner, for a comparison with OSRC preconditioning + +````julia +@target calderon_preconditioner (geo,;κ) -> begin + Γ = geo.Γ + X = raviartthomas(Γ) + Y = BEAST.buffachristiansen(Γ) + + T = Maxwell3D.singlelayer(wavenumber=κ) + N = NCross() + + Tyy = assemble(T,Y,Y); + Nxy = Matrix(assemble(N,X,Y)); + iNxy = inv(Nxy); + P = iNxy' * Tyy * iNxy + return (;P=P) +end + +@target solution_cald_EFIE (formulation_EFIE, calderon_preconditioner; residual) -> begin + (;bilforms, linforms) = formulation_EFIE + (;A) = bilforms; (;bx) = linforms; + P_calderon = calderon_preconditioner.P + iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500, left_preconditioner=P_calderon) + u, ch = BEAST.solve(iT, bx) + return (;iters=ch.iters, u=u) +end +```` + +```` +0x3b488b9b8fce1697 +```` + +Set the simulation parameters and run for different mesh sizes ``h`` + +````julia +@target benchmark_OSRC (solution_EFIE, solution_OSRC_EFIE, solution_cald_EFIE, ;) -> begin + return (;iters_EFIE = solution_EFIE.iters, iters_OSRC_EFIE = solution_OSRC_EFIE.iters, iters_cald_EFIE = solution_cald_EFIE.iters) +end +@sweep benchmark_OSRC_sweep (!benchmark_OSRC,; h=[], κ=[], residual=[], Np=[], curvature=[]) -> benchmark_OSRC + +h_values = [0.3, 0.2, 0.15] +κ = 1.0*pi +residual = 1e-6 +Np = 4 +curvature = 1.0 + +df = make(benchmark_OSRC_sweep; h=h_values, κ=κ, residual=residual, Np=Np, curvature=curvature) +```` + +```@raw html +
3×8 DataFrame
Rowiters_EFIEiters_cald_EFIEhresidualcurvatureκNpiters_OSRC_EFIE
Int64?Int64?Float64?Float64?Float64?Float64?Int64?Int64?
115690.151.0e-61.03.14159448
212990.21.0e-61.03.14159447
3106100.31.0e-61.03.14159445
+``` + +Finally, the GMRES iterations for the different methods are plotted + +````julia +using Plots +plot(df.h, df.iters_EFIE, label="EFIE", l=2) +plot!(df.h, df.iters_OSRC_EFIE, label="OSRC-EFIE", l=2) +plot!(df.h, df.iters_cald_EFIE, label="Calderón-EFIE", l=2) +xlabel!("h") +ylabel!("GMRES iterations") +```` + +```@raw html + +``` + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/examples/OSRC_efie.jl b/examples/OSRC_efie.jl new file mode 100644 index 00000000..0ecf4170 --- /dev/null +++ b/examples/OSRC_efie.jl @@ -0,0 +1,128 @@ +# In this example, the preconditioning ability of the On-Surface Radiation Condition (OSRC) Magnetic-to-Electric (MtE) map is shown on the Electric Field Integral Equation (EFIE), on the sphere. +# # OSRC operators +# The MtE map ``\text{V}^{+}`` is a map of the trace of the magnetic field ``n \times h`` to the trace of the electric field ``n \times e`` on a scattering surface. +# The approximation of the local MtE surface operator as an On-Surface Radiation Condition (OSRC) operator on an arbitrary surface ``\Gamma`` is given by: +# ```math +# \begin{equation} +# \text{V}^{+} \left( n \times h \right) + n \times e = 0, \quad \text{on} \ \Gamma, +# \end{equation} +# ``` +# In the implementation a highly accurate pseudo-differential operator ``\text{V}_{\epsilon}^{+}`` is used to represent the Magnetic-to-Electric map. +# ```math +# \begin{aligned} +# \text{V}_{\epsilon}^{+} &:= - \pmb{\Theta}^{-1} \Lambda_{2, \epsilon}^{-1} \Lambda_{1, \epsilon}, \\ +# \pmb{\Theta} \left( \cdot \right) &:= \left( n \times \cdot \right), \\ +# \Lambda_{2, \epsilon} &:= I - \textbf{curl}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{curl}_{\Gamma}, \\ +# \Lambda_{1, \epsilon} &:= \left(I + \mathcal{J} \right)^{1/2}, \\ +# \mathcal{J} &:=\frac{\Delta_{\Gamma}}{\kappa_{\epsilon}^2} = \textbf{Grad}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{Div}_{\Gamma} - \textbf{curl}_{\Gamma} \frac{1}{\kappa_{\epsilon}^2} \textup{curl}_{\Gamma}. +# \end{aligned} +# ``` +# More details on this OSRC operator and its implementation are given in +# [An OSRC Preconditioner for the EFIE (Betcke et al. (2021))](https://arxiv.org/abs/2111.10761) +# The Electric-to-Magnetic (EtM) OSRC operator is also implemented, as given in +# [Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell’s equations (C. Geuzaine et al. (2014))](https://www.researchgate.net/publication/261636307_Approximate_local_magnetic-to-electric_surface_operators_for_time-harmonic_Maxwell's_equations) +# +# # OSRC EFIE preconditioning example +# We begin with loading needed packages +using BEAST, CompScienceMeshes + +using Makeitso +using LinearAlgebra + +using Plots +using PlotlyDocumenter + +# Now, assemble and solve the unpreconditioned EFIE + +BLAS.set_num_threads(8) + +@target geo (;h) -> begin + (; Γ = CompScienceMeshes.meshsphere(radius=1.0, h=h)) +end + +@target formulation_EFIE (geo,; κ) -> begin + X = raviartthomas(geo.Γ) + T = Maxwell3D.singlelayer(wavenumber=κ) + + E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ); + e = (n × E) × n; + + bx = assemble(e, X) + A = assemble(T,X,X); + return (;bilforms=(;A), linforms=(;bx)) +end + +@target solution_EFIE (formulation_EFIE,; residual) -> begin + (;bilforms, linforms) = formulation_EFIE + (;A) = bilforms; (;bx) = linforms; + iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500) + u, ch = BEAST.solve(iT, bx) + return (;iters=ch.iters, u=u) +end + +# Next, assemble the OSRC MtE operator on the primal grid and use it to precondition the EFIE + +@target OSRC_MtE_op (geo,;κ, Np, curvature) -> begin + MtE_OSRC_operator = BEAST.MtE_OSRC_op(κ, Np, pi/2, curvature) + Nd = BEAST.nedelec(geo.Γ) + MtE_map = assemble(MtE_OSRC_operator, Nd, Nd) + return (;MtE=MtE_map) +end + +@target solution_OSRC_EFIE (formulation_EFIE, OSRC_MtE_op; residual) -> begin + (;bilforms, linforms) = formulation_EFIE + (;A) = bilforms; (;bx) = linforms; + P_OSRC = OSRC_MtE_op.MtE + iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500, left_preconditioner=P_OSRC) + u, ch = BEAST.solve(iT, bx) + return (;iters=ch.iters, u=u) +end + +# Finally, assemble and use the Calderón preconditioner, for a comparison with OSRC preconditioning + +@target calderon_preconditioner (geo,;κ) -> begin + Γ = geo.Γ + X = raviartthomas(Γ) + Y = BEAST.buffachristiansen(Γ) + + T = Maxwell3D.singlelayer(wavenumber=κ) + N = NCross() + + Tyy = assemble(T,Y,Y); + Nxy = Matrix(assemble(N,X,Y)); + iNxy = inv(Nxy); + P = iNxy' * Tyy * iNxy + return (;P=P) +end + +@target solution_cald_EFIE (formulation_EFIE, calderon_preconditioner; residual) -> begin + (;bilforms, linforms) = formulation_EFIE + (;A) = bilforms; (;bx) = linforms; + P_calderon = calderon_preconditioner.P + iT = BEAST.GMRESSolver(A; restart=1_500, reltol=residual, maxiter=1_500, left_preconditioner=P_calderon) + u, ch = BEAST.solve(iT, bx) + return (;iters=ch.iters, u=u) +end + +# Set the simulation parameters and run for different mesh sizes ``h`` + +@target benchmark_OSRC (solution_EFIE, solution_OSRC_EFIE, solution_cald_EFIE, ;) -> begin + return (;iters_EFIE = solution_EFIE.iters, iters_OSRC_EFIE = solution_OSRC_EFIE.iters, iters_cald_EFIE = solution_cald_EFIE.iters) +end +@sweep benchmark_OSRC_sweep (!benchmark_OSRC,; h=[], κ=[], residual=[], Np=[], curvature=[]) -> benchmark_OSRC + +h_values = [0.3, 0.2, 0.15] +κ = 1.0*pi +residual = 1e-6 +Np = 4 +curvature = 1.0 + +df = make(benchmark_OSRC_sweep; h=h_values, κ=κ, residual=residual, Np=Np, curvature=curvature) + +# Finally, the GMRES iterations for the different methods are plotted +using Plots +plot(df.h, df.iters_EFIE, label="EFIE", l=2) +plot!(df.h, df.iters_OSRC_EFIE, label="OSRC-EFIE", l=2) +plot!(df.h, df.iters_cald_EFIE, label="Calderón-EFIE", l=2) +xlabel!("h") +ylabel!("GMRES iterations") \ No newline at end of file diff --git a/src/BEAST.jl b/src/BEAST.jl index bde48f83..b8953c3c 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -243,6 +243,7 @@ include("dyadicop.jl") include("interpolation.jl") include("operators/projectors.jl") include("operators/theta.jl") +include("operators/OSRC.jl") include("quadrature/rules/momintegrals.jl") include("quadrature/doublenumints.jl") diff --git a/src/operators/OSRC.jl b/src/operators/OSRC.jl new file mode 100644 index 00000000..13429e5b --- /dev/null +++ b/src/operators/OSRC.jl @@ -0,0 +1,253 @@ +using SparseArrays +using LinearAlgebra +using BEAST +import Polynomials + +struct MtE_OSRC_op <: Operator + wavenumber::Float64 + Np::Int + θ_p::Float64 + curvature::Float64 + aj::Vector{ComplexF64} + bj::Vector{ComplexF64} + Aj::Vector{ComplexF64} + Bj::Vector{ComplexF64} +end + +# First a rotating branch-cut rational Padé approximation of the square root function ``\sqrt{1+z^2}`` is implemented. +imag_conv = -im # The paper uses mathematical time-harmonic convention ``e^{-i \omega t}`` -> swap to match BEAST time-harmonic convention ``e^{i \omega t}``. +function MtE_OSRC_op(wavenumber::Float64, Np::Int, θ_p::Float64, curvature::Float64) + # get the real and rotated pade coefficients + aj = ComplexF64[] + bj = ComplexF64[] + Aj = ComplexF64[] + Bj = ComplexF64[] + for j in 1:Np + a = 2/(2*Np+1)*(sin(j*pi/(2*Np+1)))^2 + b = (cos(j*pi/(2*Np+1)))^2 + A = exp(-imag_conv*θ_p/2) * a / (1 + b*(exp(-imag_conv * θ_p) - 1))^2 + B = exp(-imag_conv*θ_p) * b / (1 + b*(exp(-imag_conv * θ_p) - 1)) + + push!(aj, a) + push!(bj, b) + push!(Aj, A) + push!(Bj, B) + end + return MtE_OSRC_op(wavenumber, Np, θ_p, curvature, aj, bj, Aj, Bj) +end + +function scalartype(op::MtE_OSRC_op) + T = scalartype(op.wavenumber) + return Complex{T} +end + +function get_RNp(z, OSRC) + R_Np = 1 + sum(OSRC.aj[j]*z/(1+OSRC.bj[j]*z) for j in 1:OSRC.Np) + return R_Np +end + +function get_C0(OSRC) + C0 = exp(imag_conv*OSRC.θ_p/2) * get_RNp((exp(-imag_conv*OSRC.θ_p)-1), OSRC) + return C0 +end + +function get_R0(OSRC) + C0 = get_C0(OSRC) + R0 = C0 + sum(OSRC.Aj[j]/OSRC.Bj[j] for j in 1:OSRC.Np) + return R0 +end + +function rotated_pade(z, OSRC) + return get_R0(OSRC)*I - sum(OSRC.Aj[j]*I/(OSRC.Bj[j]*(I + OSRC.Bj[j]*z)) for j in 1:OSRC.Np) +end + +function rotated_implicit_pade(get_Π, OSRC) + return get_R0(OSRC)*I - sum(OSRC.Aj[j]*I/(OSRC.Bj[j]*(get_Π(j, OSRC))) for j in 1:OSRC.Np) +end + + +# # Projection and embedding of LinearMaps +# Construct LinearMaps which efficiently slice out a relevant submatrix of a LinearMap (without constructing the actual LinearMap). + +struct SlicedLinearMap + A::LinearMap + P::LinearMap + Q::LinearMap + rows::UnitRange{Int} + cols::UnitRange{Int} +end + +function SlicedLinearMap(A::LinearMap, rows::UnitRange{Int}, cols::UnitRange{Int}) + # selected rows and columns inside the matrix + n_rows = length(rows) + n_cols = length(cols) + n, m = size(A) + + # functions used for the construction of the matrices + function slice_row(vector::AbstractVector) + return vector[rows] + end + + function slice_column(vector::AbstractVector) + y = zeros(ComplexF64, n) + y[cols] = vector + return y + end + + P = LinearMap(slice_row, n_rows, n; ismutating=false) + Q = LinearMap(slice_column, n, n_cols; ismutating=false) + return SlicedLinearMap(A, P, Q, rows, cols) +end + +# The square root operator is regularized by adding a small imaginary component ``\epsilon`` to the wavenumber: ``k_{\epsilon} = k + i \epsilon``. +function MtE_damping(op) + return 0.39*op.wavenumber^(1/3)*op.curvature^(2/3) +end + +function assemble(op::MtE_OSRC_op,X::Space,Y::Space; quadstrat=defaultquadstrat) + R_0 = get_R0(op) + + ϵ = MtE_damping(op) + κ = op.wavenumber + κ_ϵ = κ + imag_conv*ϵ + + #create auxilary basis functions + L0_int = BEAST.lagrangec0d1(X.geo) + grad_L0_int = BEAST.gradient(L0_int) + # Define the relevant function spaces + curl_X = BEAST.curl(X) + curl_Y = BEAST.curl(Y) + + N_L0 = numfunctions(L0_int) + N_X = numfunctions(X) + N_Y = numfunctions(Y) + + # Assemble the submatrices of the blockmatrix of the system + Id = BEAST.Identity(); + G = assemble(Id, X, Y; quadstrat) + N_ϵ = (1/κ_ϵ)^2 * assemble(Id, curl_X, curl_Y; quadstrat) + K_ϵ = κ_ϵ^2 * assemble(Id, L0_int, L0_int; quadstrat) + L = assemble(Id, X, grad_L0_int; quadstrat) + L_transpose = assemble(Id, grad_L0_int, Y; quadstrat) + + # construct the sparse system matrix and invert + function create_j_phi_matrix17(j) + B_j = op.Bj[j] + # blockmatrix of sparse matrices + AXY = [G-B_j*N_ϵ B_j*L + L_transpose K_ϵ] + SXY = BEAST.lu(AXY) + Sliced_SXY = SlicedLinearMap(SXY, 1:N_X, 1:N_Y) + P = Sliced_SXY.P + Q = Sliced_SXY.Q + SXY_sliced = P*SXY*Q + return SXY_sliced + end + + sum_Π_inv_matrix = sum(op.Aj[j]/op.Bj[j] * create_j_phi_matrix17(j) for j in 1:op.Np) + G_N_ϵ_inv = BEAST.lu(G - N_ϵ) + + MtE_map = - (G_N_ϵ_inv * R_0 - G_N_ϵ_inv * G * sum_Π_inv_matrix) + return MtE_map +end + +## The Electric-to-Magnetic (EtM) OSRC operator is also implemented. +# More details for this OSRC operator are given in the paper Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell’s equations (C. Geuzaine et al. (2014)). + +# get Polynomial representation of numerator and denumerator Padé series, to obtain an inversion of the fraction. + +function prod_poles(op::MtE_OSRC_op, range) + return prod(Polynomials.Polynomial([1, op.Bj[j]]) for j in range) +end + +function A_j_polynomial(op::MtE_OSRC_op, j) + range = filter(x -> x != j, 1:op.Np) + return Polynomials.Polynomial([0, op.Aj[j]]) * prod_poles(op, range) +end + +function get_numerator(op::MtE_OSRC_op) + return BEAST.get_C0(op) * prod_poles(op, 1:op.Np) + sum(A_j_polynomial(op, j) for j in 1:op.Np) +end + +function get_denumerator(op::MtE_OSRC_op) + return prod_poles(op, 1:op.Np) +end + +struct EtM_OSRC_op + wavenumber::Float64 + Np::Int + θ_p::Float64 + curvature::Float64 + r0::ComplexF64 + rj::Vector{ComplexF64} + qj::Vector{ComplexF64} +end + +# construct inverse pade rational +function EtM_OSRC_op(wavenumber::Float64, Np::Int, θ_p::Float64, curvature::Float64) + # construct MtE operator with regular rotational branch cut Padé approximation + MtE_OSRC_op = BEAST.MtE_OSRC_op(wavenumber, Np, θ_p, curvature) + pade_inv_num = get_denumerator(MtE_OSRC_op) + pade_inv_denum = get_numerator(MtE_OSRC_op) + qj = Polynomials.roots(pade_inv_denum) + + # polynmial division to get constant term r0 + r0, remainder = Polynomials.divrem(pade_inv_num, pade_inv_denum) + P_deriv = Polynomials.derivative(pade_inv_denum) + residues = [remainder(p)/P_deriv(p) for p in qj] # compute residues=poles at every point + rj = residues + return EtM_OSRC_op(wavenumber, Np, θ_p, curvature, r0, rj, qj) +end + +function rational_pade_inv(z, op::EtM_OSRC_op) + s = op.r0 + for (r, q) in zip(op.rj, op.qj) + s += r / (z - q) + end + return s +end + + +function assemble(op::EtM_OSRC_op,X::Space,Y::Space; quadstrat=defaultquadstrat) + ϵ = BEAST.MtE_damping(op) + κ = op.wavenumber + κ_ϵ = κ + imag_conv*ϵ + + # assemble M2 matrix + + curl_X = BEAST.curl(X); + curl_Y = BEAST.curl(Y); + P1 = BEAST.lagrangec0d1(X.geo); + P1_grad = BEAST.gradient(P1); + + N_X = numfunctions(X) + N_Y = numfunctions(Y) + + Id = BEAST.Identity(); + A = assemble(Id, X, Y; quadstrat); + N = (1/κ_ϵ)^2 * assemble(Id, curl_X, curl_Y; quadstrat); + K = κ_ϵ^2 * assemble(Id, P1, P1; quadstrat); + L = assemble(Id, X, P1_grad; quadstrat); + L_transpose = assemble(Id, P1_grad, Y; quadstrat) + + function get_M2(j) + Bq = - (op.qj[j]*A + N); + M2 = [Bq L + L_transpose K]; + M2_inv = BEAST.lu(M2) + SliceMap = BEAST.SlicedLinearMap(M2_inv, 1:N_X, 1:N_Y) + P = SliceMap.P + Q = SliceMap.Q + M2_inv_sliced = P*M2_inv*Q + return M2_inv_sliced + end + + A_inv = BEAST.lu(A) + EtM_map = A_inv * (op.r0*I + A * sum(op.rj[j]*get_M2(j) for j in 1:op.Np))*(A - N) + return EtM_map +end + +function scalartype(op::EtM_OSRC_op) + T = scalartype(op.wavenumber) + return Complex{T} +end \ No newline at end of file