If I did not miss something, the type MWSingleLayerField3D lacks the alpha and beta members that MWSingleLayer3D has. This is not critical as for most practical purposes it boils down to scaling with the impedance. It is, however, a mild inconsistency with MWSingleLayer3D.

I would propose to

• make MWSingleLayerField3D fully compatible with MWSingleLayer3D
• add a function that allows constructing an instance of MWSingleLayerField3D from an instance of MWSingleLayer3D (this also for the double layer)

If approved of I am happy to do this.

Consider this code

using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

Γ = meshsphere(1.0, 0.8)
X = raviartthomas(Γ)

sol = 299792458.0
dt, Nt = 3/sol,200

dirac = timebasisdelta(dt, Nt)
pulse = timebasiscxd0(dt, Nt)

duration, delay, amplitude = 150.0/sol, 300.0/sol, 1.0
direction, polarisation = z, x
gaussian = creategaussian(duration, delay, amplitude)
gaussian_delayed = creategaussian(duration, delay+dt, amplitude)

E = planewave(polarisation, direction, gaussian, sol)
E_delayed = planewave(polarisation, direction, gaussian_delayed, sol)
Ederiv = planewave(polarisation, direction, derive(gaussian), sol)

Vreference = assemble(Ederiv, X ⊗ pulse)
Vtheory = assemble(E, X ⊗ dirac) - assemble(E_delayed, X ⊗ dirac)

There might be 2 issues with the pulse (timebasiscxd0) testing:

• it seems that there is a delay of 1 time step as the first column of Vreference is all 0
• Considering the following derivation i think that Vreference should be equal to Vtheory but instead there is a dt scaling in this code and effectively Vreference[:,2:end]/dt equals Vtheory[:,1:end-1] ([:,2:end] and [:,1:end-1] because of the delay of 1 time step described in the first point)

Hi! In the readme it says that the efie.jl is for a perfectly conducting sphere. Could you please tell me where this enters the code and whether it can be changed?

Thank you very much!

Hello, I would like to try your package but I cannot installed it on Julia 1.0.3:

(v1.0) pkg> add https://github.com/krcools/BEAST.jl
   Cloning git-repo `https://github.com/krcools/BEAST.jl`
  Updating git-repo `https://github.com/krcools/BEAST.jl`
 Resolving package versions...
ERROR: Unsatisfiable requirements detected for package StaticArrays [90137ffa]:
 StaticArrays [90137ffa] log:
 ├─possible versions are: [0.0.1-0.0.11, 0.1.0-0.1.5, 0.2.0-0.2.1, 0.3.0-0.3.1, 0.4.0, 0.5.0-0.5.1, 0.6.0-0.6.6, 0.7.0-0.7.2, 0.8.0-0.8.3, 0.9.0-0.9.2, 0.10.0, 0.10.2] or uninstalled
 ├─restricted to versions 0.8.3-0.9 by BEAST [bb4162c7], leaving only versions [0.8.3, 0.9.0-0.9.2]
 │ └─BEAST [bb4162c7] log:
 │   ├─possible versions are: 1.0.0 or uninstalled
 │   └─BEAST [bb4162c7] is fixed to version 1.0.0+
 └─restricted to versions 0.10.2 by an explicit requirement — no versions left

Perhaps a loose of version requirements would work? It probably involves some compatible test. Thanks!

Certain mesh discretization sizes for certain geometries e.g simply chaging tdefie.jl - Line 8 to
**Δx = 1/10**
Γ = meshrectangle(1.0,1.0,Δx)
results in

LoadError: inconsistent contour construction!
in WiltonInts84

Using Δx = 1/8 or Δx = 1/12 for example doesn't return any error.

@varform η*T[k,j] == f[ntrace(k,γ₁)]

returns error

MethodError: objects of type Mesh{3,2,Float64} are not callable

Stacktrace:
 [1] top-level scope at none:0

Now positions need to be supplied every time a RTSpace or LagSpace is created. This propagates to taking the divergence and trace operations on bases.

In addition positions do not make sense for global basis functions.

Solution probably is to introduce many more types for the bases.

In addition the meaning of the trace and divergence and how this interacts with matrix indices needs to be redesigned.

Running a statement like below
RT = rt_ports(Γ,[γ₁,γ₂])
f=ScalarTrace(p->1.0)
assemble(f,ntrace(RT,γ₁))

results in errors like:

objects of type BEAST.ScalarTrace are not callable
MethodError: no method matching integrand(::BEAST.ScalarTrace,...

That’s probably because Line 143 and 144 of helmholtop.jl have been commented out.
Line 143
Line 144

If a quadrature point happens to be on a vertex or edge of a triangle of the refinement, finite Cauchy integrals are split into multiple infinite contributions, breaking the assemble routine. This only happens in practice if the trial mesh is a refinement of the test mesh.

Trying to run the helmhotz3d_neumann.jl example, after fixing the path and using issue, it stops at

julia> x1 = solve(eq1)
ERROR: MethodError: no method matching (BlockArrays.PseudoBlockArray{ComplexF64, N, R, BS} where {N, R<:AbstractArray{ComplexF64, N}, BS<:Tuple{Vararg{AbstractUnitRange{Int64}, N}}})(::Vector{ComplexF64}, ::Vector{Int64})
Closest candidates are:
  (BlockArrays.PseudoBlockArray{T1, N, R, BS} where {N, R<:AbstractArray{T1, N}, BS<:Tuple{Vararg{AbstractUnitRange{Int64}, N}}})(::AbstractArray{T2, N}) where {T1, T2, N} at C:\Users\pi96doc\.julia\packages\BlockArrays\sRJ26\src\pseudo_blockarray.jl:104
  (BlockArrays.PseudoBlockArray{T, N, R, BS} where {N, R<:AbstractArray{T, N}, BS<:Tuple{Vararg{AbstractUnitRange{Int64}, N}}})(::UndefInitializer, ::AbstractVector{Int64}...) where {T, N} at C:\Users\pi96doc\.julia\packages\BlockArrays\sRJ26\src\pseudo_blockarray.jl:76
Stacktrace:
 [1] assemble(lform::BEAST.Variational.LinForm, test_space_dict::Dict{Any, Any})
   @ BEAST ~\Documents\Programming\Julia\BEAST.jl\src\solvers\solver.jl:129
 [2] solve(eq::BEAST.DiscreteEquation)
   @ BEAST ~\Documents\Programming\Julia\BEAST.jl\src\solvers\lusolver.jl:19
 [3] top-level scope
   @ REPL[16]:1

Any ideas how to fix this?

Some issue finding length of a StaticArray when using weld in CompScienceMeshes after migrating to BEAST v1.0.0. Stacktrace:

I am currently working on problems that require me to assemble BilForms with a large number of terms. And even for small numbers of unknowns, the assembly took way longer than I had expected.

This seems to be because in the assembly, the LinearMaps are added by a '+', which results in a tuple, which if I understand correctly then causes long compile times. @sbadrian fixed this once in commit 989f285, but was reverted in commit f8784ed.

Now, I have tried to fix this in PR #134, and for my case it works, and tests are still running, but I feel like I have butchered the code a bit in the process, as there now are if-statements checking whether it's a SpaceTimeBasis also before the return.

If you have any idea how to improve what I did there or how to do it differently, I would be really grateful!

I tried running an example with the current GitHub version and hit a few issues:
I noted that the path for the examples together with reading some repository-contained data is a bit cumbersome to use in practice as it requires a cd() to the examples folder.
How about starting the examples like this:

d = dirname(pathof(BEAST)) * "/../examples"
Γ = readmesh(joinpath(d,"sphere2.in"))

It turns out, Julia in Windows can deal with the "/" as path separators, so this is tested in Windows ...
in addition you may want to add

using LinearAlgebra

as the norm function is not part of julia Base.

For mixed discretizations involving both primal and dual basis functions, blockassembler function delivers results deviating from assembler function.

PR #89 updates blockassembler function to replicate the assembler function way of handling these discretizations.

Testing the double layer operator K with nxRWG and the BuffaChristiansen as the trial basis results in NaN values for certain interactions.
e.g

RT = raviartthomas(Γ); BC = buffachristiansen(Γ)
K = MWDoubleLayer3D(im*κ)
M  = assemble(K,RT,BC)

Tried with both open and closed structures

Currently, modified (real-valued gamma) and general Helmholtz kernels (complex gamma) are supported.

However, there are dedicated fast methods for laplace kernels (gamma = 0.0), which currently could not be distinguished from a modified Helmholtz kernel (since 0.0 is a real-valued number).

I was wondering what could be a way out. For example, integer kernels could be considered as static (but that leads to contradictions, if gamma = 1).

I would suggest to reserve gamma = nothing for the case that it is static. In fact, this allows to efficiently dispatch on some operations avoiding a multiplication with 0.

In the long run, if more general kernels are implemented, this is probably the soundest way to check whether an extension exists (in which case it is not nothing).

PR #103 implements this feature request.

Note: The last commit of PR #103 contains a bugfix that allows to use complex DirectProductSpaces (currently, they result in excessive compile times since in the background tuples are used). I can put this in another pull request if desired.

Recently, we realised that the sign of the DoubleLayerNear Helmholtz3D operator doesn't match with the literature, e.g. Sauter & Schwab - in PR #91 this is fixed and occurances of the operator are adapted accordingly.

new position vector in RTBasis type RTBasis
causes this freq domain code to fail
RT = rt_ports(Γ,[γ₁,γ₂])
T = MWSingleLayer3D(im*κ)
Z = assemble(T,RT,RT)
because rtports has no implementation for populating the position vector yet. Which is accessed by
subset <- assemble!

I will try and include this later, although I wonder how the position value will be averaged for the "global" component of rtports as in Line 54

Using the single layer time domain operator appropriately with the speed of light,
and a suitable time step for a simple example like in tdefie.jl where:
T = MWSingleLayerTDIO(c, -1/c, -c, 2, 0) or T = MWSingleLayerTDIO(c)
Δt, Nt = 3e-10, 200
Results in current coefficients of the order 1e-47, growing with time, and eventually NaN values. Same applies to structures of different geometry.

I was wondering if there is a reason why the type can be provided as a keyword argument. I think in principle is suffices to use a promote_type on scalartype(op) and eltype(coeffs).

If there is nothing withstanding, I can change the behavior.

Hi, I am Jürgen Furhmann and work on VoronoiFVM.
I am organizing a minisymposium at JuliaCon 2020, see https://github.com/users/j-fu/projects/1
Would you be interested ?

The blockassembler mechanism for building system matrices seems to deliver results deviating from the standard assemble interface, when mixed discretizations are used, for example, BC testing and RWG expansion functions (i.e., I did not yet test the Helmholtz/Laplace case). This was tested with single-threaded assembly.

My tests indicate that the standard assemble routines delivers results that can be considered as the “correct” ones.

`using BEAST
using Test
using CompScienceMeshes
using LinearAlgebra

function hassemble(operator::BEAST.AbstractOperator,
test_functions,
trial_functions)

blkasm = BEAST.blockassembler(operator, test_functions, trial_functions)

function assembler(Z, tdata, sdata)
    store(v,m,n) = (Z[m,n] += v)
    blkasm(tdata,sdata,store)
end

mat = zeros(scalartype(operator), 
            numfunctions(test_functions), 
            numfunctions(trial_functions))

assembler(mat, 1:numfunctions(test_functions), 1:numfunctions(trial_functions))
return mat

end

c = 3e8
μ = 4π1e-7
ε = 1/(μc^2)
f = 1e8
λ = c/f
k = 2π/λ
ω = k*c
η = sqrt(μ/ε)

a = 1
Γ = CompScienceMeshes.meshcuboid(a,a,a,0.2)

𝓣 = Maxwell3D.singlelayer(wavenumber=k)
𝓚 = Maxwell3D.doublelayer(wavenumber=k)

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

println("Number of RWG functions: ", numfunctions(X))

T_blockassembler = hassemble(𝓣, X, X)
T_standardassembler = assemble(𝓣, X, X)

@test norm(T_blockassembler - T_standardassembler)/norm(T_standardassembler) ≈ 0.0 atol=1e-14

T_bc_blockassembler = hassemble(𝓣, Y, Y)
T_bc_standardassembler = assemble(𝓣, Y, Y)

@test norm(T_bc_blockassembler - T_bc_standardassembler)/norm(T_bc_standardassembler) ≈ 0.0 atol=1e-14

K_mix_blockassembler = hassemble(𝓚,Y,X)
K_mix_standardassembler = assemble(𝓚,Y,X)

T_mix_blockassembler = hassemble(𝓣, Y, X)
T_mix_standardassembler = assemble(𝓣, Y, X)

@test norm(K_mix_blockassembler - K_mix_standardassembler)/norm(K_mix_standardassembler) ≈ 0.0 atol=1e-14
@test norm(T_mix_blockassembler - T_mix_standardassembler)/norm(T_mix_standardassembler) ≈ 0.0 atol=1e-14

`

I am sorry to put the request here but I can not get in touch with you in any other way. I am applying a grant ask me a proof that I have participated in previous projects. Can you help me provide a letter to confirm that I have involved in the following project?
2015.8 - 2018.1 EPSRC project: Boundary Element Methods for Next-Gen Devices in TeraHertz Technology. (University of Nottingham)
If you need any information, please feel free to contact me. email address: [email protected].

Zhaowei Liu

Hello,
I tried to run the helmholtz3d_dirichlet.jl example and I received the following error:

ERROR: LoadError: MethodError: no method matching +(::StaticArrays.SVector{3, ComplexF64}, ::ComplexF64)
For element-wise addition, use broadcasting with dot syntax: array .+ scalar
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
  +(::BEAST.Polynomial, ::Number) at /home/rsclinux/.julia/packages/BEAST/Y3rc9/src/utils/polynomial.jl:97
  +(::ChainRulesCore.Tangent{P, T} where T, ::P) where P at /home/rsclinux/.julia/packages/ChainRulesCore/Y1Mee/src/tangent_arithmetic.jl:145
  ...
Stacktrace:
 [1] facecurrents(coeffs::Vector{ComplexF64}, basis::BEAST.LagrangeBasis{0, -1, Mesh{3, 3, Float64}, Float64, 1, StaticArrays.SVector{3, Float64}})
   @ BEAST ~/.julia/packages/BEAST/Y3rc9/src/postproc.jl:37
 [2] top-level scope
   @ ~/Documents/Playing-Around/Beast/main.jl:31
in expression starting at /home/rsclinux/Documents/Playing-Around/Beast/main.jl:31

It this something related to the package?

For this analysis, I am considering a single pair of triangles and as function spaces Raviart-Thomas (X) and Buffa-Christansen (Y).

For the X-Y (test-ansatz) combination, I observe extremely slow convergences. As quadstrat, I have employed DoubleNumWiltonSauterQStrat. Strangely, only vertex-Sauter-Schwab hits are detected in the default quadrule-routine.

The code used for testing is

using StaticArrays
using CompScienceMeshes
using BEAST
using SauterSchwabQuadrature
using LinearAlgebra

CM= CompScienceMeshes

"""
    generate_refpair()

Generates a reference pair of triangles to test the implementation of Nedelec elements.
"""
function generate_refpair(;angle=180)
    vertices = [
        SVector(0.0, 1.0, 0.0),
        SVector(-0.5, 0.0, 0.0),
        SVector(+0.5, 0.0, 0.0),
        SVector(0.0, 1.0*cos(angle/180*π), 1.0*sin(angle/180*π))
    ]
    
    triangles = [
        SVector(1, 2, 3),
        SVector(2, 4, 3)
    ]
    return Mesh(vertices, triangles)
end
Γ = generate_refpair(angle=45)

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

Kop = Maxwell3D.doublelayer(wavenumber=0.0)

for NPoint in (NPoints = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40])
    K = Matrix(
            assemble(
                Kop,
                X,
                Y,
                quadstrat=BEAST.DoubleNumWiltonSauterQStrat(1,1,1,1,10,10,10,NPoint)
            )
        )
    @show NPoint
    @show K[1,1]
end

It results in something that is, at best, only slowly converging (looking NPoint = 20, 30, 40) :

NPoint = 1
K[1, 1] = 0.3133105682715872
NPoint = 2
K[1, 1] = 0.08489650498161906
NPoint = 3
K[1, 1] = 0.3022434482571144
NPoint = 4
K[1, 1] = 0.13065516181679618
NPoint = 5
K[1, 1] = 0.28213562023208744
NPoint = 6
K[1, 1] = 0.14725373021626786
NPoint = 7
K[1, 1] = 0.268205868252778
NPoint = 8
K[1, 1] = 0.15485095040328895
NPoint = 9
K[1, 1] = 0.2578434780988771
NPoint = 10
K[1, 1] = 0.15995038734757924
NPoint = 15
K[1, 1] = 0.23834860397099786
NPoint = 20
K[1, 1] = 0.17026104335668432
NPoint = 30
K[1, 1] = 0.17381558608911996
NPoint = 40
K[1, 1] = 0.1756178103848572

The X-X case results, as expected, in both edge and vertex Sauter-Schwab cases. The evaluation of

for NPoint in (NPoints = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40])
    K = Matrix(
            assemble(
                Kop,
                X,
                X,
                quadstrat=BEAST.DoubleNumWiltonSauterQStrat(1,1,1,1,10,NPoint,NPoint,10)
            )
        )
    @show NPoint
    @show K[1,1]
    push!(Ks, K[1,1])
end

results in

NPoint = 1
K[1, 1] = 9.839311909382144e-18
NPoint = 2
K[1, 1] = 1.8716953419283273e-18
NPoint = 3
K[1, 1] = 1.4136733354495269e-18
NPoint = 4
K[1, 1] = 1.514786870765978e-18
NPoint = 5
K[1, 1] = 4.296681196888337e-20
NPoint = 6
K[1, 1] = 2.0987591709721173e-18
NPoint = 7
K[1, 1] = 3.358793160516325e-19
NPoint = 8
K[1, 1] = 4.273795823654268e-19
NPoint = 9
K[1, 1] = -5.456270038354023e-19
NPoint = 10
K[1, 1] = 1.0845503997686032e-18
NPoint = 15
K[1, 1] = 2.853070397214177e-19
NPoint = 20
K[1, 1] = 8.91422364348399e-19
NPoint = 30
K[1, 1] = 3.8998795267911214e-19
NPoint = 40
K[1, 1] = -1.4806170598824172e-19

I did not verify it, but probably the result should be exactly zero, and then this might be due to limitations of finite precision.
Edit: Double checked it in Mathematica and for this case, it should be exactly zero, so the X-X appears to be fine.

In test/test_hh3d_nearfield, the spherical mesh generated by gmsh is different on WIndows and linux for the specific parameters passed to CompScienceMeshes.meshsphere. This can result in error bounds set in a Linux environment to be too tight for Windows.

It would be better to load mesh files from an asset subdirectory of test than to generate the mesh on every run. Not only is this more robust, it also results in fewer non-test related CPU time.

Currently, there is only a plane-wave excitation available for the Helmholtz3D kernel. For electrostatic simulations a uniform electric field (i.e., a linear potential) would be interesting. In addition, both for the static and dynamic regime, a monopole source is useful (and in future, also a dipole source could be added).

While implementing, I had the impression that there was no Dirichlet-Trace available. #87 adds this and the aforementioned excitations.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!