Merge branch 'master' of github.com:krcools/BEAST.jl

src/BEAST.jl

@@ -224,6 +224,7 @@ include("helmholtz3d/nitsche.jl")
 include("helmholtz3d/hh3dnear.jl")
 include("helmholtz3d/hh3dfar.jl")
 include("helmholtz3d/hh3d_sauterschwabqr.jl")
+include("helmholtz3d/helmholtz3d.jl")
 
 #suport for Volume Integral equation
 include("volumeintegral/vie.jl")

src/decoupled/dpops.jl

@@ -1,4 +1,4 @@
-struct CurlSingleLayerDP3D{T,U} <: MaxwellOperator3D
+struct CurlSingleLayerDP3D{T,U} <: MaxwellOperator3D{U,T}
     gamma::T
     alpha::U
 end

src/helmholtz3d/helmholtz3d.jl

@@ -0,0 +1,120 @@
+module Helmholtz3D
+
+    using ..BEAST
+    const Mod = BEAST
+    using LinearAlgebra
+
+    function singlelayer(;
+        alpha=nothing,
+        gamma=nothing,
+        wavenumber=nothing
+    )
+
+        alpha, gamma = Mod.operator_parameter_handler(alpha, gamma, wavenumber)
+
+        return Mod.HH3DSingleLayerFDBIO(alpha,gamma)
+    end
+
+    function doublelayer(;
+        alpha=nothing,
+        gamma=nothing,
+        wavenumber=nothing
+    )
+
+        alpha, gamma = Mod.operator_parameter_handler(alpha, gamma, wavenumber)
+
+        return Mod.HH3DDoubleLayerFDBIO(alpha, gamma)
+    end
+
+    function doublelayer_transposed(;
+        alpha=nothing,
+        gamma=nothing,
+        wavenumber=nothing
+    )
+
+        alpha, gamma = Mod.operator_parameter_handler(alpha, gamma, wavenumber)
+
+        return Mod.HH3DDoubleLayerTransposedFDBIO(alpha, gamma)
+    end
+
+    function hypersingular(;
+        alpha=nothing,
+        beta=nothing,
+        gamma=nothing,
+        wavenumber=nothing
+    )
+
+        gamma, wavenumber = Mod.gamma_wavenumber_handler(gamma, wavenumber)
+
+        if alpha === nothing
+            if gamma !== nothing
+                alpha = gamma^2
+            else
+                alpha = 0.0 # In the long run, this should probably be rather 'nothing'
+            end
+        end
+
+        if beta === nothing
+            if gamma !== nothing
+                beta = one(gamma)
+            else
+                beta = one(alpha)
+            end
+        end
+
+        return Mod.HH3DHyperSingularFDBIO(alpha, beta, gamma)
+    end
+
+    function planewave(;
+            direction=error("direction is a required argument"),
+            gamma=nothing,
+            wavenumber=nothing,
+            amplitude=one(eltype(direction)))
+
+        gamma, wavenumber = Mod.gamma_wavenumber_handler(gamma, wavenumber)
+
+        # Note: Unlike for the operators, there seems little benefit in
+        # explicitly declaring a Laplace-Type excitation.
+
+        gamma === nothing && (gamma = zero(amplitude))
+
+        return Mod.HH3DPlaneWave(direction, amplitude, gamma)
+    end
+
+    function linearpotential(; direction=SVector(1, 0, 0), amplitude=1.0)
+        return Mod.HH3DLinearPotential(direction ./ norm(direction), amplitude)
+    end
+
+    function grad_linearpotential(; direction=SVector(0.0, 0.0, 0.0), amplitude=1.0)
+        return Mod.gradHH3DLinearPotential(direction, amplitude)
+    end
+
+    function monopole(;
+        position=SVector(0.0, 0.0, 0.0),
+        gamma=nothing,
+        wavenumber=nothing,
+        amplitude=1.0
+    )
+
+        gamma, wavenumber = Mod.gamma_wavenumber_handler(gamma, wavenumber)
+        gamma === nothing && (gamma = zero(amplitude))
+
+        return Mod.HH3DMonopole(position, gamma, amplitude)
+    end
+
+    function grad_monopole(;
+        position=SVector(0.0, 0.0, 0.0),
+        gamma=nothing,
+        wavenumber=nothing,
+        amplitude=1.0
+    )
+
+        gamma, wavenumber = Mod.gamma_wavenumber_handler(gamma, wavenumber)
+        gamma === nothing && (gamma = zero(amplitude))
+
+        return Mod.gradHH3DMonopole(position, gamma, amplitude)
+    end
+
+end
+
+export Helmholtz3D

src/helmholtz3d/hh3d_sauterschwabqr.jl

@@ -14,7 +14,7 @@ function pulled_back_integrand(op::HH3DSingleLayerFDBIO,
         j = jacobian(x) * jacobian(y)
 
         α = op.alpha
-        γ = op.gamma
+        γ = gamma(op)
         R = norm(cartesian(x)-cartesian(y))
         G = exp(-γ*R)/(4*π*R)
 
@@ -45,7 +45,7 @@ function pulled_back_integrand(op::HH3DHyperSingularFDBIO,
 
         α = op.alpha
         β = op.beta
-        γ = op.gamma
+        γ = gamma(op)
         R = norm(cartesian(x)-cartesian(y))
         G = exp(-γ*R)/(4*π*R)
 
@@ -85,7 +85,7 @@ function pulled_back_integrand(op::HH3DDoubleLayerFDBIO,
         j = jacobian(x) * jacobian(y)
 
         α = op.alpha
-        γ = op.gamma
+        γ = gamma(op)
 
         r = cartesian(x) - cartesian(y)
         R = norm(r)
@@ -118,7 +118,7 @@ function pulled_back_integrand(op::HH3DDoubleLayerTransposedFDBIO,
         j = jacobian(x) * jacobian(y)
 
         α = op.alpha
-        γ = op.gamma
+        γ = gamma(op)
 
         r = cartesian(x) - cartesian(y)
         R = norm(r)

src/helmholtz3d/hh3dexc.jl

@@ -1,57 +1,42 @@
-
-
-struct HH3DPlaneWave{T,P}
+struct HH3DPlaneWave{P,K,T}
     direction::P
-    wavenumber::T
+    gamma::K
     amplitude::T
 end
 
-function HH3DPlaneWave(direction, wavenumber; amplitude=1)
-    w, a = promote(wavenumber, amplitude)
-    HH3DPlaneWave(direction, w, a)
-end
-
 function (f::HH3DPlaneWave)(r)
     d = f.direction
-    k = f.wavenumber
+    γ = f.gamma
     a = f.amplitude
-    a * exp(-im*k*dot(d,r))
+    return a * exp(-γ*dot(d,r))
 end
 
-scalartype(f::HH3DPlaneWave{T}) where {T} = complex(T)
+scalartype(f::HH3DPlaneWave{P,K,T}) where {P,K,T} = promote_type(eltype(P), K, T)
 
 """
     HH3DLinearPotential
 
 A potential that linearly increases in `direction` with scaling coefficient `amplitude`.
 Its negative gradient will be a uniform vector field pointing in the opposite direction.
 """
-struct HH3DLinearPotential{T,P}
+struct HH3DLinearPotential{P,T}
     direction::P
     amplitude::T
 end
 
-function HH3DLinearPotential(; direction=SVector(1,0,0), amplitude=1.0)
-    HH3DLinearPotential(direction ./ norm(direction), amplitude)
-end
+scalartype(f::HH3DLinearPotential{P,T}) where {P,T} = promote_type(eltype(P), T)
 
 function (f::HH3DLinearPotential)(r)
     d = f.direction
     a = f.amplitude
     return a * dot(d, r)
 end
 
-scalartype(f::HH3DLinearPotential{T}) where {T} = complex(T)
-
 struct gradHH3DLinearPotential{T,P}
     direction::P
     amplitude::T
 end
 
-function gradHH3DLinearPotential(;direction=SVector(0.0,0.0,0.0), amplitude=1.0)
-    gradHH3DLinearPotential(direction, amplitude)
-end
-
 function (f::gradHH3DLinearPotential)(r)
     d = f.direction
     a = f.amplitude
@@ -73,58 +58,48 @@ dot(::NormalVector, m::gradHH3DLinearPotential) = NormalDerivative(HH3DLinearPot
 
 Potential of a monopole-type point source (e.g., of an electric charge)
 """
-struct HH3DMonopole{T,P}
+struct HH3DMonopole{P,K,T}
     position::P
-    wavenumber::T
+    gamma::K
     amplitude::T
 end
 
-scalartype(x::HH3DMonopole{T}) where {T} = complex(T)
-
-function HH3DMonopole(;position=SVector(0.0,0.0,0.0), wavenumber=0.0, amplitude=1.0)
-    w, a = promote(wavenumber, amplitude)
-    HH3DMonopole(position, w, a)
-end
+scalartype(x::HH3DMonopole{P,K,T}) where {P,K,T} = promote_type(eltype(P), K, T)
 
 function (f::HH3DMonopole)(r)
-    k = f.wavenumber
+    γ = f.gamma
     p = f.position
     a = f.amplitude
 
-    return  a*exp(-im * k * norm(r - p)) / (norm(r - p))
+    return  a*exp(-γ * norm(r - p)) / (norm(r - p))
 end
 
-struct gradHH3DMonopole{T,P}
+struct gradHH3DMonopole{P,K,T}
     position::P
-    wavenumber::T
+    gamma::K
     amplitude::T
 end
 
-scalartype(x::gradHH3DMonopole{T}) where {T} = complex(T)
-
-function gradHH3DMonopole(;position=SVector(0.0,0.0,0.0), wavenumber=0.0, amplitude=1.0)
-    w, a = promote(wavenumber, amplitude)
-    gradHH3DMonopole(position, w, a)
-end
+scalartype(x::gradHH3DMonopole{P,K,T}) where {P,K,T} = promote_type(eltype(P), K, T)
 
 function (f::gradHH3DMonopole)(r)
     a = f.amplitude
-    k = f.wavenumber
+    γ = f.gamma
     p = f.position
     vecR = r - p
     R = norm(vecR)
 
-    return -a * vecR * exp(-im * k * R) / R^2 * (im * k + 1 / R)
+    return -a * vecR * exp(-γ * R) / R^2 * (γ + 1 / R)
 end
 
 function grad(m::HH3DMonopole)
-    return gradHH3DMonopole(m.position, m.wavenumber, m.amplitude)
+    return gradHH3DMonopole(m.position, m.gamma, m.amplitude)
 end
 
-*(a::Number, m::HH3DMonopole) = HH3DMonopole(m.position, m.wavenumber, a * m.amplitude)
-*(a::Number, m::gradHH3DMonopole) = gradHH3DMonopole(m.position, m.wavenumber, a * m.amplitude)
+*(a::Number, m::HH3DMonopole) = HH3DMonopole(m.position, m.gamma, a * m.amplitude)
+*(a::Number, m::gradHH3DMonopole) = gradHH3DMonopole(m.position, m.gamma, a * m.amplitude)
 
-dot(::NormalVector, m::gradHH3DMonopole) = NormalDerivative(HH3DMonopole(m.position, m.wavenumber, m.amplitude))
+dot(::NormalVector, m::gradHH3DMonopole) = NormalDerivative(HH3DMonopole(m.position, m.gamma, m.amplitude))
 
 mutable struct DirichletTrace{T,F} <: Functional
     field::F
@@ -162,11 +137,11 @@ end
 
 function (f::NormalDerivative{T,F})(manipoint) where {T,F<:HH3DPlaneWave}
     d = f.field.direction
-    k = f.field.wavenumber
+    γ = f.field.gamma
     a = f.field.amplitude
     n = normal(manipoint)
     r = cartesian(manipoint)
-    -im*k*a * dot(d,n) * exp(-im*k*dot(d,r))
+    return -γ*a * dot(d,n) * exp(-γ*dot(d,r))
 end
 
 function (f::NormalDerivative{T,F})(manipoint) where {T,F<:HH3DLinearPotential}