diff --git a/src/Nystrom.jl b/src/Nystrom.jl index defb9e1..365a012 100644 --- a/src/Nystrom.jl +++ b/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 diff --git a/src/dim.jl b/src/dim.jl index 8f23497..03ad6ff 100644 --- a/src/dim.jl +++ b/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) diff --git a/src/integraloperator.jl b/src/integraloperator.jl index ab19e48..8638a72 100644 --- a/src/integraloperator.jl +++ b/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) diff --git a/src/kernels.jl b/src/kernels.jl index fc03266..10fa69c 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -88,21 +88,72 @@ end HyperSingularKernel{T}(op) where {T} = HyperSingularKernel{T,typeof(op)}(op) HyperSingularKernel(op) = HyperSingularKernel{default_kernel_eltype(op)}(op) +""" + struct CombinedFieldKernel{T,Op} <: AbstractKernel{T} + +Given an operator `Op`, construct its free-space combined field kernel. This +corresponds to a linear combination of `S::SingleLayerKernel` and `D::DoubleLayerKernel` +of the form `α*D - β*S`. +""" +struct CombinedFieldKernel{T,Op,Tm<:Number} <: AbstractKernel{T} + pde::Op + α::Tm + β::Tm +end +CombinedFieldKernel{T}(op,α::Tm,β::Tm) where {T,Tm} = CombinedFieldKernel{T,typeof(op),Tm}(op,α,β) +CombinedFieldKernel(op,α,β) = CombinedFieldKernel{default_kernel_eltype(op)}(op,promote(α,β)...) + +""" + struct AdjointCombinedFieldKernel{T,Op} <: AbstractKernel{T} + +Given an operator `Op`, construct its free-space adjoint combined field kernel. This +corresponds to a linear combination of `K::AdjointDoubleLayerKernel` and `H::HyperSingularKernel` +of the form `α*H - β*K`. +""" +struct AdjointCombinedFieldKernel{T,Op,Tm<:Number} <: AbstractKernel{T} + pde::Op + α::Tm + β::Tm +end +AdjointCombinedFieldKernel{T}(op,α::Tm,β::Tm) where {T,Tm} = AdjointCombinedFieldKernel{T,typeof(op),Tm}(op,α,β) +AdjointCombinedFieldKernel(op,α,β) = AdjointCombinedFieldKernel{default_kernel_eltype(op)}(op,promote(α,β)...) + # a trait for the kernel type struct SingleLayer end struct DoubleLayer end struct AdjointDoubleLayer end struct HyperSingular end +struct CombinedField end +struct AdjointCombinedField end -kernel_type(::SingleLayerKernel) = SingleLayer() -kernel_type(::DoubleLayerKernel) = DoubleLayer() -kernel_type(::AdjointDoubleLayerKernel) = AdjointDoubleLayer() -kernel_type(::HyperSingularKernel) = HyperSingular() +kernel_type(::SingleLayerKernel) = SingleLayer() +kernel_type(::DoubleLayerKernel) = DoubleLayer() +kernel_type(::AdjointDoubleLayerKernel) = AdjointDoubleLayer() +kernel_type(::HyperSingularKernel) = HyperSingular() +kernel_type(::CombinedFieldKernel) = CombinedField() +kernel_type(::AdjointCombinedFieldKernel) = AdjointCombinedField() combined_field_coefficients(::SingleLayerKernel) = (0,-1) combined_field_coefficients(::DoubleLayerKernel) = (1,0) combined_field_coefficients(::AdjointDoubleLayerKernel) = (0,-1) combined_field_coefficients(::HyperSingularKernel) = (1,0) +combined_field_coefficients(r::CombinedFieldKernel) = (r.α,r.β) +combined_field_coefficients(r::AdjointCombinedFieldKernel) = (r.α,r.β) + +# definition of combined field kernels +function (CF::CombinedFieldKernel)(target,source) + SL = SingleLayerKernel(pde(CF)) + DL = DoubleLayerKernel(pde(CF)) + α,β = combined_field_coefficients(CF) + return α*DL(target,source) - β*SL(target,source) +end + +function (CF::AdjointCombinedFieldKernel)(target,source) + K = AdjointDoubleLayerKernel(pde(CF)) + H = HyperSingularKernel(pde(CF)) + α,β = combined_field_coefficients(CF) + return α*H(target,source) - β*K(target,source) +end ################################################################################ ################################# LAPLACE ###################################### diff --git a/src/potential.jl b/src/potential.jl index 7fabcea..1cf4978 100644 --- a/src/potential.jl +++ b/src/potential.jl @@ -23,3 +23,7 @@ Base.getindex(pot::IntegralPotential,σ::AbstractVector) = (x) -> pot(σ,x) SingleLayerPotential(op::AbstractPDE,surf) = IntegralPotential(SingleLayerKernel(op),surf) DoubleLayerPotential(op::AbstractPDE,surf) = IntegralPotential(DoubleLayerKernel(op),surf) +AdjointDoubleLayerPotential(op::AbstractPDE,surf) = IntegralPotential(AdjointDoubleLayerKernel(op),surf) +HyperSingularPotential(op::AbstractPDE,surf) = IntegralPotential(HyperSingularKernel(op),surf) +CombinedFieldPotential(op::AbstractPDE,surf;α,β) = IntegralPotential(CombinedFieldKernel(op,α,β),surf) +AdjointCombinedFieldPotential(op::AbstractPDE,surf;α,β) = IntegralPotential(AdjointCombinedFieldKernel(op,α,β),surf) \ No newline at end of file diff --git a/test/dim_test.jl b/test/dim_test.jl index e4cbd17..7bbb61c 100644 --- a/test/dim_test.jl +++ b/test/dim_test.jl @@ -23,6 +23,7 @@ Random.seed!(1) for pde in ops T = Nystrom.default_density_eltype(pde) c = rand(T) + α,β = rand(eltype(T),2) # combined field coefficients u = (qnode) -> SingleLayerKernel(pde)(qnode,xout)*c dudn = (qnode) -> AdjointDoubleLayerKernel(pde)(qnode,xout)*c γ₀u = Density(u,mesh) @@ -42,6 +43,9 @@ Random.seed!(1) @test norm(e0,Inf) > norm(e1,Inf) @test norm(e1,Inf) < rtol end + C = CombinedFieldOperator(pde,mesh;α,β) + Cdim = Nystrom.assemble_dim(C) + @test Nystrom.materialize(Cdim) ≈ α*Nystrom.materialize(Ddim)-β*Nystrom.materialize(Sdim) # adjoint double-layer and hypersingular K = AdjointDoubleLayerOperator(pde,mesh) Kmat = K |> Matrix @@ -55,6 +59,9 @@ Random.seed!(1) @test norm(e0,Inf) > rtol @test norm(e1,Inf) < rtol end + AC = AdjointCombinedFieldOperator(pde,mesh;α,β) + ACdim = Nystrom.assemble_dim(AC) + @test Nystrom.materialize(ACdim) ≈ α*Nystrom.materialize(Hdim)-β*Nystrom.materialize(Kdim) end end @@ -74,6 +81,7 @@ Random.seed!(1) for pde in ops T = Nystrom.default_density_eltype(pde) c = rand(T) + α,β = rand(eltype(T),2) # combined field coefficients u = (qnode) -> SingleLayerKernel(pde)(xout,qnode)*c dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xout,qnode))*c γ₀u = Density(u,mesh) @@ -93,6 +101,9 @@ Random.seed!(1) @test norm(e0,Inf) > norm(e1,Inf) @test norm(e1,Inf) < rtol end + C = CombinedFieldOperator(pde,mesh;α,β) + Cdim = Nystrom.assemble_dim(C) + @test Nystrom.materialize(Cdim) ≈ α*Nystrom.materialize(Ddim)-β*Nystrom.materialize(Sdim) # adjoint double-layer and hypersingular pde isa Maxwell && continue K = AdjointDoubleLayerOperator(pde,mesh) @@ -107,6 +118,9 @@ Random.seed!(1) @test norm(e0,Inf) > norm(e1,Inf) @test norm(e1,Inf) < rtol end + AC = AdjointCombinedFieldOperator(pde,mesh;α,β) + ACdim = Nystrom.assemble_dim(AC) + @test Nystrom.materialize(ACdim) ≈ α*Nystrom.materialize(Hdim)-β*Nystrom.materialize(Kdim) end end @@ -125,6 +139,7 @@ Random.seed!(1) for pde in ops T = Nystrom.default_density_eltype(pde) c = rand(T) + α,β = rand(eltype(T),2) # combined field coefficients u = (qnode) -> SingleLayerKernel(pde)(xin,qnode)*c dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xin,qnode))*c γ₀u = Density(u,mesh) @@ -144,6 +159,9 @@ Random.seed!(1) @test norm(e0,Inf) > norm(e1,Inf) @test norm(e1,Inf) < rtol end + C = CombinedFieldOperator(pde,mesh;α,β) + Cdim = Nystrom.assemble_dim(C) + @test Nystrom.materialize(Cdim) ≈ α*Nystrom.materialize(Ddim)-β*Nystrom.materialize(Sdim) # adjoint double-layer and hypersingular K = AdjointDoubleLayerOperator(pde,mesh) Kmat = K |> Matrix @@ -157,6 +175,9 @@ Random.seed!(1) @test norm(e0,Inf) > rtol @test norm(e1,Inf) < rtol end + AC = AdjointCombinedFieldOperator(pde,mesh;α,β) + ACdim = Nystrom.assemble_dim(AC) + @test Nystrom.materialize(ACdim) ≈ α*Nystrom.materialize(Hdim)-β*Nystrom.materialize(Kdim) end end @@ -176,6 +197,7 @@ Random.seed!(1) for pde in ops T = Nystrom.default_density_eltype(pde) c = rand(T) + α,β = rand(eltype(T),2) # combined field coefficients u = (qnode) -> SingleLayerKernel(pde)(xs,qnode)*c dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xs,qnode))*c γ₀u = Density(u,mesh) @@ -195,6 +217,9 @@ Random.seed!(1) @test norm(e0,Inf) > norm(e1,Inf) @test norm(e1,Inf) < rtol end + C = CombinedFieldOperator(pde,mesh;α,β) + Cdim = Nystrom.assemble_dim(C) + @test Nystrom.materialize(Cdim) ≈ α*Nystrom.materialize(Ddim)-β*Nystrom.materialize(Sdim) # adjoint double-layer and hypersingular pde isa Maxwell && continue K = AdjointDoubleLayerOperator(pde,mesh) @@ -209,6 +234,9 @@ Random.seed!(1) @test norm(e0,Inf) > norm(e1,Inf) @test norm(e1,Inf) < rtol end + AC = AdjointCombinedFieldOperator(pde,mesh;α,β) + ACdim = Nystrom.assemble_dim(AC) + @test Nystrom.materialize(ACdim) ≈ α*Nystrom.materialize(Hdim)-β*Nystrom.materialize(Kdim) end end end diff --git a/test/integraloperator_test.jl b/test/integraloperator_test.jl index f082cea..8f490c6 100644 --- a/test/integraloperator_test.jl +++ b/test/integraloperator_test.jl @@ -5,7 +5,7 @@ using StaticArrays @testset "Basic tests" begin pde = Helmholtz(;dim=3,k=1) - @gmsh begin + Ω,M = @gmsh begin Geometry.clear_entities!() gmsh.initialize() GmshSDK.set_meshsize(0.1) @@ -15,7 +15,7 @@ using StaticArrays gmsh.model.mesh.generate(3) Ω = GmshSDK.domain(dim=3) M = GmshSDK.meshgen(Ω,dim=3) - return Ω,M + Ω,M end Γ = boundary(Ω) mesh = NystromMesh(view(M,Γ),order=1) diff --git a/test/nystrommesh_test.jl b/test/nystrommesh_test.jl index 8fe0039..a24d19f 100644 --- a/test/nystrommesh_test.jl +++ b/test/nystrommesh_test.jl @@ -13,7 +13,7 @@ using Nystrom gmsh.model.mesh.generate(3) Ω = GmshSDK.domain(dim=3) M = GmshSDK.meshgen(Ω,dim=3) - return Ω,M + Ω,M end ∂Ω = boundary(Ω) mesh = NystromMesh(view(M,∂Ω),order=1) @@ -37,7 +37,7 @@ using Nystrom gmsh.model.mesh.generate(3) Ω = GmshSDK.domain(dim=3) M = GmshSDK.meshgen(Ω,dim=3) - return Ω,M + Ω,M end Γ = boundary(Ω) nmesh = NystromMesh(M,Γ,order=4) # NystromMesh of surface Γ @@ -66,7 +66,7 @@ using Nystrom gmsh.model.mesh.generate(2) Ω = GmshSDK.domain(dim=2) M = GmshSDK.meshgen(Ω,dim=2) - return Ω,M + Ω,M end Γ = boundary(Ω) mesh = NystromMesh(view(M,Ω),order=2)