Resolves issue #22

src/Nystrom.jl

@@ -47,12 +47,20 @@ export
     DoubleLayerKernel,
     AdjointDoubleLayerKernel,
     HyperSingularKernel,
+    CombinedFieldKernel,
+    AdjointCombinedFieldKernel,
     SingleLayerPotential,
     DoubleLayerPotential,
+    AdjointDoubleLayerPotential,
+    HyperSingularPotential,
+    CombinedFieldPotential,
+    AdjointCombinedFieldPotential,
     SingleLayerOperator,
     DoubleLayerOperator,
     AdjointDoubleLayerOperator,
     HyperSingularOperator,
+    CombinedFieldOperator,
+    AdjointCombinedFieldOperator,
     Density,
     TangentialDensity,
     # functions

src/dim.jl

@@ -4,19 +4,25 @@ function assemble_dim(iop::IntegralOperator;compress=Matrix,location=:onsurface)
     pde  = iop.kernel.pde
     T    = kernel_type(iop)
     if T === SingleLayer()
-        singlelayer_dim(pde,X,Y;compress,location)
+        return singlelayer_dim(pde,X,Y;compress,location)
     elseif T === DoubleLayer()
-        doublelayer_dim(pde,X,Y;compress,location)
+        return doublelayer_dim(pde,X,Y;compress,location)
     elseif T === AdjointDoubleLayer()
-        adjointdoublelayer_dim(pde,X,Y;compress,location)
+        return adjointdoublelayer_dim(pde,X,Y;compress,location)
     elseif T === HyperSingular()
-        hypersingular_dim(pde,X,Y;compress,location)
+        return hypersingular_dim(pde,X,Y;compress,location)
+    elseif T === CombinedField() 
+        α,β = combined_field_coefficients(iop)
+        return combinedfield_dim(pde,X,Y;α,β,compress,location)
+    elseif T === AdjointCombinedField() 
+        α,β = combined_field_coefficients(iop)
+        return adjointcombinedfield_dim(pde,X,Y;α,β,compress,location)
     else
         notimplemented()
     end
 end
 
-function single_doublelayer_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
+function _single_doublelayer_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
     msg = "unrecognized value for kw `location`: received $location.
            Valid options are `:onsurface`, `:inside` and `:outside`."
     σ = location === :onsurface ? -0.5 : location === :inside ? 0 : location === :outside ? -1 : error(msg)
@@ -38,22 +44,34 @@ function single_doublelayer_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
     # compute corrections
     @timeit_debug "compute dim correction" begin
         if pde isa Maxwell
-            δS = _singular_weights_dim_maxwell(Sop,γ₀B,γ₁B,R,dict_near)
-            δD = _singular_weights_dim_maxwell(Dop,γ₀B,γ₁B,R,dict_near)
+            δS,δD = _singular_weights_dim_maxwell(Sop,Dop,γ₀B,γ₁B,R,dict_near)
         else
-            δS = _singular_weights_dim(Sop,γ₀B,γ₁B,R,dict_near)
-            δD = _singular_weights_dim(Dop,γ₀B,γ₁B,R,dict_near)
+            δS,δD = _singular_weights_dim(Sop,Dop,γ₀B,γ₁B,R,dict_near)
         end
     end
+    return S,D,δS,δD
+end
+function single_doublelayer_dim(args...;kwargs...)
+    S,D,δS,δD = _single_doublelayer_dim(args...;kwargs...)
     # convert to DiscreteOp
-    S_discrete = DiscreteOp(S) + DiscreteOp(δS)
-    D_discrete = DiscreteOp(D) + DiscreteOp(δD)
-    return S_discrete,D_discrete
+    Sd = DiscreteOp(S) + DiscreteOp(δS)
+    Dd = DiscreteOp(D) + DiscreteOp(δD)
+    return Sd,Dd
 end
 singlelayer_dim(args...;kwargs...) = single_doublelayer_dim(args...;kwargs...)[1]
 doublelayer_dim(args...;kwargs...) = single_doublelayer_dim(args...;kwargs...)[2]
+function combinedfield_dim(pde,X,Y;α,β,compress,location)
+    Cop = CombinedFieldOperator(pde,X,Y;α,β)
+    # convert to a possibly more efficient format
+    @timeit_debug "assemble combined field dense part" begin
+        C = compress(Cop)
+    end
+    _,_,δS,δD = _single_doublelayer_dim(pde,X,Y;compress,location)
+    # convert to DiscreteOp
+    return DiscreteOp(C) + DiscreteOp(α*δD-β*δS)
+end
 
-function adjointdoublelayer_hypersingular_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
+function _adjointdoublelayer_hypersingular_dim(pde,X,Y=X;compress=Matrix,location=:onsurface)
     msg = "unrecognized value for kw `location`: received $location.
     Valid options are `:onsurface`, `:inside` and `:outside`."
     σ = location === :onsurface ? -0.5 : location === :inside ? 0 : location === :outside ? -1 : error(msg)
@@ -67,27 +85,27 @@ function adjointdoublelayer_hypersingular_dim(pde,X,Y=X;compress=Matrix,location
     basis,γ₁_basis = _basis_dim(Kop)
     γ₀B,γ₁B,R      = _auxiliary_quantities_dim(Kop,K,H,basis,γ₁_basis,σ)
     # compute corrections
-    δK        = _singular_weights_dim(Kop,γ₀B,γ₁B,R,dict_near)
-    δH        = _singular_weights_dim(Hop,γ₀B,γ₁B,R,dict_near)
+    δK,δH = _singular_weights_dim(Kop,Hop,γ₀B,γ₁B,R,dict_near)
+    return K,H,δK,δH
+end
+function adjointdoublelayer_hypersingular_dim(args...;kwargs...)
+    K,H,δK,δH = _adjointdoublelayer_hypersingular_dim(args...;kwargs...)
     # convert to DiscreteOp
-    K_discrete = DiscreteOp(K) + DiscreteOp(δK)
-    H_discrete = DiscreteOp(H) + DiscreteOp(δH)
-    return K_discrete,H_discrete
+    Kd = DiscreteOp(K) + DiscreteOp(δK)
+    Hd = DiscreteOp(H) + DiscreteOp(δH)
+    return Kd,Hd
 end
 adjointdoublelayer_dim(args...;kwargs...)  = adjointdoublelayer_hypersingular_dim(args...;kwargs...)[1]
 hypersingular_dim(args...;kwargs...)       = adjointdoublelayer_hypersingular_dim(args...;kwargs...)[2]
-
-
-function singular_weights_dim(iop::IntegralOperator,compress=Matrix)
-    X,Y,op = iop.X, iop.Y, iop.kernel.op
-    σ = X == Y ? -0.5 : 0.0
-    #
-    dict_near = near_interaction_list(dofs(X),Y;atol=0)
-    basis,γ₁_basis = _basis_dim(iop)
-    Op1, Op2       = _auxiliary_operators_dim(iop,compress)
-    γ₀B,γ₁B,R      = _auxiliary_quantities_dim(iop,Op1,Op2,basis,γ₁_basis,σ)
-    Sp = _singular_weights_dim(iop,γ₀B,γ₁B,R,dict_near)
-    return Sp # a sparse matrix
+function adjointcombinedfield_dim(pde,X,Y;α,β,compress,location)
+    Cop = AdjointCombinedFieldOperator(pde,X,Y;α,β)
+    # convert to a possibly more efficient format
+    @timeit_debug "assemble adjoint combined field dense part" begin
+        C = compress(Cop)
+    end
+    _,_,δK,δH = _adjointdoublelayer_hypersingular_dim(pde,X,Y;compress,location)
+    # convert to DiscreteOp
+    return DiscreteOp(C) + DiscreteOp(α*δH-β*δK)
 end
 
 function _auxiliary_quantities_dim(iop,Op0,Op1,basis,γ₁_basis,σ)
@@ -107,7 +125,6 @@ function _auxiliary_quantities_dim(iop,Op0,Op1,basis,γ₁_basis,σ)
     # integrate the basis over Y
     xnodes      = dofs(X)
     num_targets = length(xnodes)
-    num_sources = length(ynodes)
     R = Matrix{T}(undef,num_targets,num_basis)
     @timeit_debug "integrate dim basis" begin
         @threads for k in 1:num_basis
@@ -152,79 +169,102 @@ function _basis_dim(iop)
     return basis,γ₁_basis
 end
 
-function _singular_weights_dim(iop::IntegralOperator,γ₀B,γ₁B,R,dict_near)
-    X,Y = target_surface(iop), source_surface(iop)
-    T   = eltype(iop)
-    num_basis = size(γ₀B,2)
-    a,b = combined_field_coefficients(iop)
-    # we now have the residue R. For the correction we need the coefficients.
+function _singular_weights_dim(op1::IntegralOperator,op2::IntegralOperator,γ₀B,γ₁B,R,dict_near)
+    # initialize vectors for the sparse matrix, then dispatch to type-stable
+    # method for each element type
+    @assert (kernel_type(op1) isa SingleLayer && kernel_type(op2) isa DoubleLayer) ||
+            (kernel_type(op1) isa AdjointDoubleLayer && kernel_type(op2) isa HyperSingular)
+    Y  = source_surface(op1)
+    T  = eltype(op1)
+    sizeop = size(op1)
     Is = Int[]
     Js = Int[]
-    Vs = T[]
+    Ss = T[]  # for single layer / adjoint double layer
+    Ds = T[]  # for double layer / hypersingular
     for (E,list_near) in dict_near
-        el2qnodes = elt2dof(Y,E)
-        num_qnodes, num_els   = size(el2qnodes)
-        M                     = Matrix{T}(undef,2*num_qnodes,num_basis)
-        @assert length(list_near) == num_els
-        for n in 1:num_els
-            j_glob                = @view el2qnodes[:,n]
-            M[1:num_qnodes,:]     = @view γ₀B[j_glob,:]
-            M[num_qnodes+1:end,:] = @view γ₁B[j_glob,:]
-            # distinguish scalar and vectorial case
+        _singular_weights_dim!(Is,Js,Ss,Ds,Y,γ₀B,γ₁B,R,E,list_near)
+    end
+    # for single layer / adjoint double layer
+    _,b = combined_field_coefficients(op1)
+    Sp = sparse(Is,Js,b*Ss,sizeop...) 
+    # for double layer / hypersingular
+    a,_ = combined_field_coefficients(op2)
+    Dp = sparse(Is,Js,a*Ds,sizeop...) 
+    return Sp, Dp
+end
+
+@noinline function _singular_weights_dim!(Is,Js,Ss,Ds,Y,γ₀B,γ₁B,R,E,list_near)
+    T = eltype(Ss)
+    num_basis = size(γ₀B,2)
+    el2qnodes = elt2dof(Y,E)
+    num_qnodes, num_els   = size(el2qnodes)
+    M                     = Matrix{T}(undef,2*num_qnodes,num_basis)
+    @assert length(list_near) == num_els
+    for n in 1:num_els
+        j_glob                = @view el2qnodes[:,n]
+        M[1:num_qnodes,:]     = @view γ₀B[j_glob,:]
+        M[num_qnodes+1:end,:] = @view γ₁B[j_glob,:]
+        # distinguish scalar and vectorial case
+        if T <: Number
+            F                     = qr(M)
+        elseif T <: SMatrix
+            M_mat = blockmatrix_to_matrix(M)
+            F                     = qr!(M_mat)
+        else
+            error("unknown element type T=$T")
+        end
+        for (i,_) in list_near[n]
             if T <: Number
-                F                     = qr(M)
+                tmp = ((R[i:i,:])/F.R)*adjoint(F.Q)
             elseif T <: SMatrix
-                M_mat = blockmatrix_to_matrix(M)
-                F                     = qr!(M_mat)
+                tmp_scalar  = (blockmatrix_to_matrix(R[i:i,:])/F.R)*adjoint(F.Q)
+                tmp  = matrix_to_blockmatrix(tmp_scalar,T)
             else
                 error("unknown element type T=$T")
             end
-            for (i,_) in list_near[n]
-                if T <: Number
-                    tmp = ((R[i:i,:])/F.R)*adjoint(F.Q)
-                elseif T <: SMatrix
-                    tmp_scalar  = (blockmatrix_to_matrix(R[i:i,:])/F.R)*adjoint(F.Q)
-                    tmp  = matrix_to_blockmatrix(tmp_scalar,T)
-                else
-                    error("unknown element type T=$T")
-                end
-                w    = axpby!(a,view(tmp,1:num_qnodes),b,view(tmp,(num_qnodes+1):(2*num_qnodes)))
-                append!(Is,fill(i,num_qnodes))
-                append!(Js,j_glob)
-                append!(Vs,w)
-            end
+            append!(Is,fill(i,num_qnodes))
+            append!(Js,j_glob)
+            append!(Ss,view(tmp,(num_qnodes+1):(2*num_qnodes)))
+            append!(Ds,view(tmp,1:num_qnodes))
         end
     end
-    Sp = sparse(Is,Js,Vs,size(iop)...)
-    return Sp
+    return Is,Js,Ss,Ds
 end
 
-function _singular_weights_dim_maxwell(iop::IntegralOperator,γ₀B,γ₁B,R,dict_near)
+function _singular_weights_dim_maxwell(Sop::IntegralOperator,Dop::IntegralOperator,γ₀B,γ₁B,R,dict_near)
     # initialize vectors for the sparse matrix, then dispatch to type-stable
     # method for each element type
-    T   = eltype(iop)
+    Y = source_surface(Sop)
+    T  = eltype(Sop)
+    sizeop = size(Sop)
     Is = Int[]
     Js = Int[]
-    Vs = T[]
+    Ss = T[]  # for single layer
+    Ds = T[]  # for double layer
     for (E,list_near) in dict_near
-        _singular_weights_dim_maxwell!(Is,Js,Vs,iop,γ₀B,γ₁B,R,E,list_near)
+        _singular_weights_dim_maxwell!(Is,Js,Ss,Ds,Y,γ₀B,γ₁B,R,E,list_near)
     end
-    Sp = sparse(Is,Js,Vs,size(iop)...)
-    return Sp
+    # for single layer
+    _,b = combined_field_coefficients(Sop)
+    Sp = sparse(Is,Js,b*Ss,sizeop...) 
+    # for double layer
+    a,_ = combined_field_coefficients(Dop)
+    Dp = sparse(Is,Js,a*Ds,sizeop...) 
+    return Sp, Dp
 end
 
-@noinline function _singular_weights_dim_maxwell!(Is,Js,Vs,iop,γ₀B,γ₁B,R,E,list_near)
-    X,Y    = target_surface(iop), source_surface(iop)
+@noinline function _singular_weights_dim_maxwell!(Is,Js,Ss,Ds,Y,γ₀B,γ₁B,R,E,list_near)
     qnodes = dofs(Y)
-    T   = eltype(iop)
-    a,b = combined_field_coefficients(iop)
     el2qnodes = elt2dof(Y,E)
     num_qnodes, num_els   = size(el2qnodes)
     num_basis = size(γ₀B,2)
+    T = eltype(Ss)
     T2 = SMatrix{2,3,ComplexF64,6}
     T3 = SMatrix{3,2,ComplexF64,6}
     M      = Matrix{T2}(undef,2*num_qnodes,num_basis)
     M_mat  = Matrix{ComplexF64}(undef,2*2*num_qnodes,3*num_basis)
+    ws     = Vector{T}(undef,num_qnodes)  # for single layer
+    wd     = Vector{T}(undef,num_qnodes)  # for double layer
     @assert length(list_near) == num_els
     for n in 1:num_els
         j_glob                = @view el2qnodes[:,n]
@@ -242,19 +282,18 @@ end
         for (i,_) in list_near[n]
             tmp_scalar  = (blockmatrix_to_matrix(R[i:i,:])/F.R)*adjoint(F.Q)
             tmp         = matrix_to_blockmatrix(tmp_scalar,T3)
-            tmp         = axpby!(a,view(tmp,1:num_qnodes),b,view(tmp,(num_qnodes+1):(2*num_qnodes)))
-            w           = Vector{T}(undef,num_qnodes)
             for k in 1:num_qnodes
-                # J,_ = jacobian(qnodes[j_glob[k]]) |> qr
                 J    = jacobian(qnodes[j_glob[k]])
-                w[k] = tmp[k]*transpose(J)
+                wd[k] = tmp[k]*transpose(J)
+                ws[k] = tmp[k+num_qnodes]*transpose(J)
             end
             append!(Is,fill(i,num_qnodes))
             append!(Js,j_glob)
-            append!(Vs,w)
+            append!(Ss,ws)
+            append!(Ds,wd)
         end
     end
-    return Is,Js,Vs
+    return Is,Js,Ss,Ds
 end
 
 function _source_gen(iop::IntegralOperator,kfactor=5)

src/integraloperator.jl

@@ -44,6 +44,8 @@ SingleLayerOperator(op::AbstractPDE,X,Y=X)        = IntegralOperator(SingleLayer
 DoubleLayerOperator(op::AbstractPDE,X,Y=X)        = IntegralOperator(DoubleLayerKernel(op), X, Y)
 AdjointDoubleLayerOperator(op::AbstractPDE,X,Y=X) = IntegralOperator(AdjointDoubleLayerKernel(op), X, Y)
 HyperSingularOperator(op::AbstractPDE,X,Y=X)      = IntegralOperator(HyperSingularKernel(op), X, Y)
+CombinedFieldOperator(op::AbstractPDE,X,Y=X;α,β)  = IntegralOperator(CombinedFieldKernel(op,α,β), X, Y)
+AdjointCombinedFieldOperator(op::AbstractPDE,X,Y=X;α,β)  = IntegralOperator(AdjointCombinedFieldKernel(op,α,β), X, Y)
 
 ambient_dimension(iop::IntegralOperator) = ambient_dimension(iop.kernel)