krcools / beast.jl Goto Github PK
View Code? Open in Web Editor NEWBoundary Element Analysis and Simulation Toolkit
License: Other
Boundary Element Analysis and Simulation Toolkit
License: Other
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!
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 ?
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!
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.
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?
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
@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
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.
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.
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.
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.
We are currently implementing the general p-order case for the Lagrange elements. For p=2 there is a (partial) implementation that implies a certain ordering of the local functions:
function (f::LagrangeRefSpace{T,2,3})(t) where T
u,v,w, = barycentric(t)
j = jacobian(t)
p = t.patch
#curl=(p[3]-p[2])/j),
# (value=v, curl=(p[1]-p[3])/j),
# (value=w, curl=(p[2]-p[1])/j)
σ = sign(dot(normal(t), cross(p[1]-p[3],p[2]-p[3])))
SVector(
(value=u*(2*u-1), curl=σ*(p[3]-p[2])*(4u-1)/j),
(value=v*(2*v-1), curl=σ*(p[1]-p[3])*(4v-1)/j),
(value=w*(2*w-1), curl=σ*(p[2]-p[1])*(4w-1)/j),
(value=4*v*w, curl=4*σ*(w*(p[1]-p[3])+v*(p[2]-p[1]))/j),
(value=4*w*u, curl=4*σ*(w*(p[3]-p[2])+u*(p[2]-p[1]))/j),
(value=4*u*v, curl=4*σ*(u*(p[1]-p[3])+v*(p[3]-p[2]))/j),
)
end
It seems that the first three functions are on the vertices and then come the functions on the edges.
For the general p-case, it is more convenient to use the multiindex notation (i,j,k) (see Figure 4.3 in [1] P. P. Silvester and R. L. Ferrari, Finite elements for electrical engineers. Cambridge University Press, 1996). Silvester use different conversion strategy multi-index -> single-index, basically, by starting from index (3,0,0), then reducing the first index to 2 and setting the second index to the maximum (2,1,0) etc.
This conversion strategy would require to change the currently used ordering. I am happy to keep the existing one, but it is not clear to me what the general pth order case is in the current ordering system.
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!
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!
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.
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.
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:
Vreference
is all 0Vreference
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)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
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?
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 matchingintegrand(::BEAST.ScalarTrace,...
That’s probably because Line 143 and 144 of helmholtop.jl have been commented out.
Line 143
Line 144
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.
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
`
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.
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
If approved of I am happy to do this.
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
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.