From 48765b614b72018b81224128e95e1e3a732157ed Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Mon, 17 Jul 2023 16:34:39 +0200
Subject: [PATCH 01/24] adding quadrule for assemble(nedelec space,double
layer, raviarthomas)
---
src/bases/local/ndlocal.jl | 23 +++++++++++++++++++++++
src/maxwell/mwops.jl | 32 ++++++++++++++++++++++++++++++++
src/maxwell/nxdbllayer.jl | 5 +++--
src/quadrature/quadrule.jl | 24 ++++++++++++++++++++++++
4 files changed, 82 insertions(+), 2 deletions(-)
diff --git a/src/bases/local/ndlocal.jl b/src/bases/local/ndlocal.jl
index 9f0c196f..f4308306 100644
--- a/src/bases/local/ndlocal.jl
+++ b/src/bases/local/ndlocal.jl
@@ -67,3 +67,26 @@ function restrict(ϕ::NDRefSpace{T}, dom1, dom2) where T
return Q
end
+
+const _vert_perms_nd = [
+ (1,2,3),
+ (2,3,1),
+ (3,1,2),
+ (2,1,3),
+ (1,3,2),
+ (3,2,1),
+]
+const _dof_perms_nd = [
+ (1,2,3),
+ (3,1,2),
+ (2,3,1),
+ (2,1,3),
+ (1,3,2),
+ (3,2,1),
+]
+
+function dof_permutation(::NDRefSpace, vert_permutation)
+ i = findfirst(==(tuple(vert_permutation...)), _vert_perms_nd)
+ @assert i != nothing
+ return _dof_perms_nd[i]
+end
\ No newline at end of file
diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl
index 0cb87dec..21b2b061 100644
--- a/src/maxwell/mwops.jl
+++ b/src/maxwell/mwops.jl
@@ -136,6 +136,38 @@ function quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j
qd.bpoints[1,j],)
end
+function quadrule(op::MaxwellOperator3D, g::NDRefSpace, f::RTRefSpace, i, τ, j, σ, qd,
+ qs::DoubleNumWiltonSauterQStrat)
+
+T = eltype(eltype(τ.vertices))
+hits = 0
+dtol = 1.0e3 * eps(T)
+dmin2 = floatmax(T)
+for t in τ.vertices
+ for s in σ.vertices
+ d2 = LinearAlgebra.norm_sqr(t-s)
+ dmin2 = min(dmin2, d2)
+ hits += (d2 < dtol)
+ end
+end
+
+hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
+hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
+hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
+
+h2 = volume(σ)
+xtol2 = 0.2 * 0.2
+k2 = abs2(op.gamma)
+max(dmin2*k2, dmin2/16h2) < xtol2 && return WiltonSERule(
+ qd.tpoints[2,i],
+ DoubleQuadRule(
+ qd.tpoints[2,i],
+ qd.bpoints[2,j],),)
+return DoubleQuadRule(
+ qd.tpoints[1,i],
+ qd.bpoints[1,j],)
+end
+
function quadrule(op::MaxwellOperator3D, g::BDMRefSpace, f::BDMRefSpace, i, τ, j, σ, qd,
qs::DoubleNumWiltonSauterQStrat)
diff --git a/src/maxwell/nxdbllayer.jl b/src/maxwell/nxdbllayer.jl
index 8772e354..e8fab643 100644
--- a/src/maxwell/nxdbllayer.jl
+++ b/src/maxwell/nxdbllayer.jl
@@ -101,7 +101,7 @@ function (igd::Integrand{<:DoubleLayerRotatedMW3D})(x,y,f,g)
# x = neighborhood(igd.test_chart,u)
# y = neighborhood(igd.trial_chart,v)
- j = jacobian(x) * jacobian(y)
+ # j = jacobian(x) * jacobian(y)
nx = normal(x)
r = cartesian(x) - cartesian(y)
@@ -117,7 +117,8 @@ function (igd::Integrand{<:DoubleLayerRotatedMW3D})(x,y,f,g)
fvalue = getvalue(f)
gvalue = getvalue(g)
- jKg = cross.(Ref(K), j*gvalue)
+ #jKg = cross.(Ref(K), j*gvalue)
+ jKg = cross.(Ref(K), gvalue)
jnxKg = cross.(Ref(nx), jKg)
return _krondot(fvalue, jnxKg)
end
diff --git a/src/quadrature/quadrule.jl b/src/quadrature/quadrule.jl
index b4eda1b7..12cd7de1 100644
--- a/src/quadrature/quadrule.jl
+++ b/src/quadrature/quadrule.jl
@@ -20,6 +20,30 @@ function quadrule(op::IntegralOperator, g::RTRefSpace, f::RTRefSpace, i, τ, j,
hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
+ return DoubleQuadRule(
+ qd.tpoints[1,i],
+ qd.bpoints[1,j],)
+end
+
+function quadrule(op::IntegralOperator, g::RTRefSpace, f::NDRefSpace, i, τ, j, σ, qd,
+ qs::DoubleNumSauterQstrat)
+
+ T = eltype(eltype(τ.vertices))
+ hits = 0
+ dtol = 1.0e3 * eps(T)
+ dmin2 = floatmax(T)
+ for t in τ.vertices
+ for s in σ.vertices
+ d2 = LinearAlgebra.norm_sqr(t-s)
+ dmin2 = min(dmin2, d2)
+ hits += (d2 < dtol)
+ end
+ end
+
+ hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
+ hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
+ hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
+
return DoubleQuadRule(
qd.tpoints[1,i],
qd.bpoints[1,j],)
From d508af7629e6f6889994ceb945a568842276e6c9 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Mon, 31 Jul 2023 18:22:47 +0200
Subject: [PATCH 02/24] update maxwell operator to any linear basis on
triangles
---
src/maxwell/mwops.jl | 5 ++++-
1 file changed, 4 insertions(+), 1 deletion(-)
diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl
index 21b2b061..a92d2b1e 100644
--- a/src/maxwell/mwops.jl
+++ b/src/maxwell/mwops.jl
@@ -104,7 +104,9 @@ end
regularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DReg(op.gamma)
singularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DSng(op.gamma)
-function quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd,
+const LinearRefSpaceTriangle = Union{RTRefSpace, NDRefSpace, BDMRefSpace, NCrossBDMRefSpace}
+
+function quadrule(op::MaxwellOperator3D, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, i, τ, j, σ, qd,
qs::DoubleNumWiltonSauterQStrat)
T = eltype(eltype(τ.vertices))
@@ -195,6 +197,7 @@ function quadrule(op::MaxwellOperator3D, g::BDMRefSpace, f::BDMRefSpace, i, τ,
end
+
function qrib(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd)
dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
From 9f65b9a8c149d6a4216e6ec61ad9c9fc9c15b1c7 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Mon, 20 Nov 2023 14:14:17 +0100
Subject: [PATCH 03/24] integration of acustic singlelayer with all analytic
strategy
---
Project.toml | 1 +
examples/tdacusticsinglelayer.jl | 51 ++++++++++
src/BEAST.jl | 1 +
src/maxwell/mwops.jl | 29 +-----
src/maxwell/timedomain/acustictdops.jl | 125 +++++++++++++++++++++++++
src/quadrature/quadstrats.jl | 3 +
src/timedomain/tdintegralop.jl | 6 ++
7 files changed, 189 insertions(+), 27 deletions(-)
create mode 100644 examples/tdacusticsinglelayer.jl
create mode 100644 src/maxwell/timedomain/acustictdops.jl
diff --git a/Project.toml b/Project.toml
index bbe15231..e2df26f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,6 +28,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SparseMatrixDicts = "5cb6c4b0-9b79-11e8-24c9-f9621d252589"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TimeDomainBEMInt = "0303f91b-5811-4998-82cd-e0687f55e8f9"
WiltonInts84 = "a3e2863e-c0ee-5ff6-a523-307a4cdc8724"
[compat]
diff --git a/examples/tdacusticsinglelayer.jl b/examples/tdacusticsinglelayer.jl
new file mode 100644
index 00000000..5d1a3f71
--- /dev/null
+++ b/examples/tdacusticsinglelayer.jl
@@ -0,0 +1,51 @@
+using CompScienceMeshes, BEAST, LinearAlgebra
+Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
+Γ = meshsphere(radius=1.0, h=0.5)
+
+X = lagrangecxd0(Γ)
+numfunctions(X)
+Δt = 0.05
+Nt = 20
+T = timebasisshiftedlagrange(Δt, Nt, 0)
+U = timebasisdelta(Δt, Nt)
+
+V = X ⊗ T
+W = X ⊗ U
+
+duration = 2 * 20 * Δt
+delay = 1.5 * duration
+amplitude = 1.0
+gaussian = creategaussian(duration, delay, amplitude)
+direction, polarisation = ẑ, x̂
+E = planewave(polarisation, direction, derive(gaussian), 1.0)
+
+@hilbertspace j
+@hilbertspace j′
+
+SL = TDAcustic3D.acusticsinglelayer(speedofsound=1.0, numdiffs=1)
+# BEAST.@defaultquadstrat (SL, W, V) BEAST.OuterNumInnerAnalyticQStrat(7)
+
+tdacusticsl = @discretise SL[j′,j] == -1.0E[j′] j∈V j′∈W
+xacusticsl = solve(tdacusticsl)
+#corregere da qui
+import Plots
+Plots.plot(xefie[1,:])
+
+import Plotly
+fcr, geo = facecurrents(xefie[:,125], X)
+Plotly.plot(patch(geo, norm.(fcr)))
+
+
+
+
+
+
+Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
+ω = collect(ω0 .+ (0:Nt-1)*Δω)
+_, i1 = findmin(abs.(ω.-1.0))
+
+ω1 = ω[i1]
+ue = Xefie[:,i1] / fouriertransform(gaussian)(ω1)
+
+fcr, geo = facecurrents(ue, X)
+Plotly.plot(patch(geo, norm.(fcr)))
\ No newline at end of file
diff --git a/src/BEAST.jl b/src/BEAST.jl
index 17e685d8..66b442fb 100644
--- a/src/BEAST.jl
+++ b/src/BEAST.jl
@@ -236,6 +236,7 @@ include("helmholtz3d/timedomain/tdhh3dops.jl")
include("helmholtz3d/timedomain/tdhh3dexc.jl")
include("helmholtz3d/timedomain/tdhh3dpp.jl")
+include("maxwell/timedomain/acustictdops.jl")#support for acusticsinglelayer
include("maxwell/timedomain/mwtdops.jl")
include("maxwell/timedomain/mwtdexc.jl")
include("maxwell/timedomain/tdfarfield.jl")
diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl
index a92d2b1e..866af4f6 100644
--- a/src/maxwell/mwops.jl
+++ b/src/maxwell/mwops.jl
@@ -106,7 +106,7 @@ singularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DSng(op.gamma)
const LinearRefSpaceTriangle = Union{RTRefSpace, NDRefSpace, BDMRefSpace, NCrossBDMRefSpace}
-function quadrule(op::MaxwellOperator3D, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, i, τ, j, σ, qd,
+function quadrule(op::MaxwellOperator3D, g::RTRefSpace, f::RTRefSpace, i, τ, j, σ, qd,
qs::DoubleNumWiltonSauterQStrat)
T = eltype(eltype(τ.vertices))
@@ -138,7 +138,7 @@ function quadrule(op::MaxwellOperator3D, g::LinearRefSpaceTriangle, f::LinearRef
qd.bpoints[1,j],)
end
-function quadrule(op::MaxwellOperator3D, g::NDRefSpace, f::RTRefSpace, i, τ, j, σ, qd,
+function quadrule(op::MaxwellOperator3D, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, i, τ, j, σ, qd,
qs::DoubleNumWiltonSauterQStrat)
T = eltype(eltype(τ.vertices))
@@ -170,31 +170,6 @@ return DoubleQuadRule(
qd.bpoints[1,j],)
end
-function quadrule(op::MaxwellOperator3D, g::BDMRefSpace, f::BDMRefSpace, i, τ, j, σ, qd,
- qs::DoubleNumWiltonSauterQStrat)
-
- hits = 0
- dtol = 1.0e3 * eps(eltype(eltype(τ.vertices)))
- dmin2 = floatmax(eltype(eltype(τ.vertices)))
- for t in τ.vertices
- for s in σ.vertices
- d2 = LinearAlgebra.norm_sqr(t-s)
- dmin2 = min(dmin2, d2)
- hits += (d2 < dtol)
- end
- end
-
- hits == 3 && return SauterSchwabQuadrature.CommonFace(qd.gausslegendre[3])
- hits == 2 && return SauterSchwabQuadrature.CommonEdge(qd.gausslegendre[2])
- hits == 1 && return SauterSchwabQuadrature.CommonVertex(qd.gausslegendre[1])
-
- h2 = volume(σ)
- xtol2 = 0.2 * 0.2
- k2 = abs2(op.gamma)
- return DoubleQuadRule(
- qd.tpoints[1,i],
- qd.bpoints[1,j],)
-end
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
new file mode 100644
index 00000000..c1f4100c
--- /dev/null
+++ b/src/maxwell/timedomain/acustictdops.jl
@@ -0,0 +1,125 @@
+mutable struct AcusticSingleLayerTDIO{T} <: RetardedPotential{T}
+ "speed of sound in medium"
+ speed_of_light::T #for compatibility with the rest of the package I left light instead of sound
+ "weight"
+ weight::T
+ "number of temporal differentiations"
+ diffs::Int
+end
+
+function Base.:*(a::Number, op::AcusticSingleLayerTDIO)
+ @info "scalar product a * op (acusticSL)"
+ AcusticSingleLayerTDIO(
+ op.speed_of_light,
+ a * op.weight,
+ op.diffs)
+end
+
+AcusticSingleLayerTDIO(;speedofsound) = AcusticSingleLayerTDIO(speedofsound, one(speedofsound), 0)
+
+module TDAcustic3D
+import ...BEAST
+
+function acusticsinglelayer(;speedofsound, numdiffs)
+ @assert numdiffs >= 0
+ #controllare con Kristof che tipo di schema temporale viene utilizzato e come si con integrali rispetto al tempo extra
+ #return BEAST.integrate(BEAST.MWSingleLayerTDIO(speedoflight,-1/speedoflight,-speedoflight,2,0))
+ return BEAST.AcusticSingleLayerTDIO(speedofsound,one(speedofsound),numdiffs)
+end
+
+end #of the module
+
+export TDAcustic3D
+
+defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1)
+
+
+
+
+
+function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
+ testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
+ return 0.0
+end
+
+quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
+ p, testel, q, trialel, r, timeel, qd, ::AllAnalyticalQStrat) = ZuccottiStrat(1.0)
+
+function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
+ testels, trialels, timeels, quadstrat::OuterNumInnerAnalyticQStrat)
+
+ dmax = numfunctions(timerefs)-1
+ bn = binomial.((0:dmax),(0:dmax)')
+
+ V = eltype(testels[1].vertices)
+ ws = WiltonInts84.workspace(V)
+ # quadpoints(testrefs, testels, (3,)), bn, ws
+ quadpoints(testrefs, testels, (quadstrat.outer_rule,)), bn, ws
+end
+
+
+quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
+ p, testel, q, trialel, r, timeel, qd, ::OuterNumInnerAnalyticQStrat) = WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])
+
+function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::ZuccottiStrat) where T
+ if ι[1]<0
+ t1=0.0
+ else
+ t1=ι[1]
+ end
+ t2=ι[2]
+
+
+ @assert t2 > t1
+
+ sos = op.speed_of_light
+
+ a1index,a2index,a3index,b1index,b2index,b3index=0,0,0,0,0,0
+ Tt = eltype(eltype(τ.vertices))
+ hits = 0
+ dtol = 1.0e3 * eps(Tt)
+ dmin2 = floatmax(Tt)
+ for t in 1:3
+ for s in 1:3
+ d2 = LinearAlgebra.norm_sqr(τ[t]-σ[s])
+ dmin2 = min(dmin2, d2)
+ if (d2 < dtol)
+ hits+=1
+ if hits==1
+ a1index =t
+ b1index = s
+ elseif hits==2
+ a2index=t
+ b2index=s
+ end
+ end
+ end
+ end
+ #print(τ[1]," ",τ[2]," ",τ[3]," ",σ[1]," ",σ[2]," ",σ[3]," ",t1," ",t2," ",)
+ if hits==3
+ z[1,1,1]+=TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2)
+ elseif hits==2
+
+ if mod1(a1index +1,3)==a2index
+ a3index=mod1(a1index-1,3)
+ else
+ a3index=mod1(a1index+1,3)
+ end
+
+ if mod1(b1index+1,3)==b2index
+ b3index=mod1(a1index-1,3)
+ else
+ b3index=mod1(a1index+1,3)
+ end
+
+ z[1,1,1]+=TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[a3index],t1,t2)
+ elseif hits==1
+ a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
+ b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
+
+ z[1,1,1]+=TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2)
+ else
+ z[1,1,1]+=inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2)
+ end
+ #for the moment sos=1 but I will correct this
+end
\ No newline at end of file
diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl
index 40d4bddb..1c583f07 100644
--- a/src/quadrature/quadstrats.jl
+++ b/src/quadrature/quadstrats.jl
@@ -38,6 +38,9 @@ struct OuterNumInnerAnalyticQStrat{R}
outer_rule::R
end
+struct AllAnalyticalQStrat{R}
+ rule::R
+end
defaultquadstrat(op, tfs, bfs) = defaultquadstrat(op, refspace(tfs), refspace(bfs))
macro defaultquadstrat(dop, body)
diff --git a/src/timedomain/tdintegralop.jl b/src/timedomain/tdintegralop.jl
index 77a9c2e7..316c274f 100644
--- a/src/timedomain/tdintegralop.jl
+++ b/src/timedomain/tdintegralop.jl
@@ -1,4 +1,5 @@
using WiltonInts84
+using TimeDomainBEMInt
abstract type RetardedPotential{T} <: Operator end
Base.eltype(::RetardedPotential{T}) where {T} = T
@@ -289,6 +290,9 @@ struct WiltonInts84Strat{T,V,W}
workspace::W
end
+struct ZuccottiStrat{T}
+ weight::T #Allanalyticalformula
+end
function momintegrals!(z, op, g, f, T, τ, σ, ι, qr::WiltonInts84Strat)
@@ -300,3 +304,5 @@ function momintegrals!(z, op, g, f, T, τ, σ, ι, qr::WiltonInts84Strat)
end
end
+
+
From a75d2f918912ac6d0c347b1d9d66296dc9366965 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Mon, 20 Nov 2023 14:52:40 +0100
Subject: [PATCH 04/24] update acustisl
---
src/maxwell/timedomain/acustictdops.jl | 15 ++++++---------
src/timedomain/tdintegralop.jl | 2 +-
2 files changed, 7 insertions(+), 10 deletions(-)
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
index c1f4100c..9bbbcca0 100644
--- a/src/maxwell/timedomain/acustictdops.jl
+++ b/src/maxwell/timedomain/acustictdops.jl
@@ -39,11 +39,11 @@ defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1)
function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
- return 0.0
+ return nothing
end
quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- p, testel, q, trialel, r, timeel, qd, ::AllAnalyticalQStrat) = ZuccottiStrat(1.0)
+ p, testel, q, trialel, r, timeel, qd, ::AllAnalyticalQStrat) = ZuccottiRule(1.0)
function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
testels, trialels, timeels, quadstrat::OuterNumInnerAnalyticQStrat)
@@ -62,12 +62,9 @@ quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
p, testel, q, trialel, r, timeel, qd, ::OuterNumInnerAnalyticQStrat) = WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])
function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::ZuccottiStrat) where T
- if ι[1]<0
- t1=0.0
- else
- t1=ι[1]
- end
- t2=ι[2]
+
+ t1=ι[1]
+ t2=ι[2]
@assert t2 > t1
@@ -99,7 +96,7 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
if hits==3
z[1,1,1]+=TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2)
elseif hits==2
-
+ #pay attention with double layer and index permutation
if mod1(a1index +1,3)==a2index
a3index=mod1(a1index-1,3)
else
diff --git a/src/timedomain/tdintegralop.jl b/src/timedomain/tdintegralop.jl
index 316c274f..3adc0d54 100644
--- a/src/timedomain/tdintegralop.jl
+++ b/src/timedomain/tdintegralop.jl
@@ -290,7 +290,7 @@ struct WiltonInts84Strat{T,V,W}
workspace::W
end
-struct ZuccottiStrat{T}
+struct ZuccottiRule{T}
weight::T #Allanalyticalformula
end
From 19c13ddb611924037d116436ded9f1b38ebe43e5 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Thu, 11 Jan 2024 13:56:09 +0100
Subject: [PATCH 05/24] update in acoustic single layer methods
---
examples/disabled/tdhh3d_neumann.jl | 20 ++-
examples/tdacusticsinglelayer.jl | 60 +++++---
src/helmholtz3d/timedomain/tdhh3dexc.jl | 1 +
src/maxwell/timedomain/acustictdops.jl | 195 +++++++++++++++++++++---
src/timedomain/tdintegralop.jl | 6 +
test/runtests.jl | 3 +-
6 files changed, 231 insertions(+), 54 deletions(-)
diff --git a/examples/disabled/tdhh3d_neumann.jl b/examples/disabled/tdhh3d_neumann.jl
index 0c7f6508..5db0c60b 100644
--- a/examples/disabled/tdhh3d_neumann.jl
+++ b/examples/disabled/tdhh3d_neumann.jl
@@ -2,8 +2,8 @@ using CompScienceMeshes
using BEAST
using LinearAlgebra
-# G = meshsphere(1.0, 0.25)
G = meshsphere(1.0, 0.30)
+#G = CompScienceMeshes.meshmobius(h=0.1)
c = 1.0
S = BEAST.HH3DSingleLayerTDBIO(c)
@@ -20,9 +20,9 @@ h = BEAST.gradient(e)
# X = lagrangecxd0(G)
X = lagrangecxd0(G)
-Y = duallagrangec0d1(G)
-# Y = lagrangecxd0(G)
-
+#Y = duallagrangec0d1(G)
+Y = lagrangecxd0(G)
+numfunctions(X)
# X = duallagrangecxd0(G, boundary(G))
# Y = lagrangec0d1(G)
@@ -31,7 +31,7 @@ Y = duallagrangec0d1(G)
Δt, Nt = 0.16, 300
# T = timebasisc0d1(Δt, Nt)
P = timebasiscxd0(Δt, Nt)
-H = timebasisc0d1(Δt, Nt)
+#H = timebasisc0d1(Δt, Nt)
δ = timebasisdelta(Δt, Nt)
# assemble the right hand side
@@ -43,9 +43,17 @@ Zd = Z0d + (-0.5)*Z1d
u = marchonintime(inv(Zd[:,:,1]), Zd, bd, Nt)
bs = assemble(e, X ⊗ δ)
-Zs = assemble(S, X ⊗ δ, X ⊗ P, Val{:bandedstorage})
+Zs = assemble(S, X ⊗ δ, X ⊗ P)
v = marchonintime(inv(Zs[:,:,1]), Zs, -bs, Nt)
+
+tdacusticsl = @discretise S[j′,j] == -1.0e[j′] j∈ (X ⊗ P) j′∈ (X ⊗ δ)
+xacusticsl = solve(tdacusticsl)
+
+
+import Plots
+Plots.plot(xacusticsl[1000,2900:3000])
+
U, Δω, ω0 = fouriertransform(u, Δt, 0.0, 2)
V, Δω, ω0 = fouriertransform(v, Δt, 0.0, 2)
diff --git a/examples/tdacusticsinglelayer.jl b/examples/tdacusticsinglelayer.jl
index 5d1a3f71..4f049814 100644
--- a/examples/tdacusticsinglelayer.jl
+++ b/examples/tdacusticsinglelayer.jl
@@ -1,23 +1,20 @@
using CompScienceMeshes, BEAST, LinearAlgebra
-Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
-Γ = meshsphere(radius=1.0, h=0.5)
+#Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
+Γ = meshsphere(radius=1.0, h=0.30)
X = lagrangecxd0(Γ)
-numfunctions(X)
-Δt = 0.05
-Nt = 20
+
+Δt = 0.16
+Nt = 300
T = timebasisshiftedlagrange(Δt, Nt, 0)
U = timebasisdelta(Δt, Nt)
V = X ⊗ T
W = X ⊗ U
-duration = 2 * 20 * Δt
-delay = 1.5 * duration
-amplitude = 1.0
-gaussian = creategaussian(duration, delay, amplitude)
-direction, polarisation = ẑ, x̂
-E = planewave(polarisation, direction, derive(gaussian), 1.0)
+width, delay, scaling = 16.0, 24.0, 1.0
+gaussian = creategaussian(width, delay, scaling)
+e = BEAST.planewave(point(0,0,1), 1.0, gaussian)
@hilbertspace j
@hilbertspace j′
@@ -25,20 +22,34 @@ E = planewave(polarisation, direction, derive(gaussian), 1.0)
SL = TDAcustic3D.acusticsinglelayer(speedofsound=1.0, numdiffs=1)
# BEAST.@defaultquadstrat (SL, W, V) BEAST.OuterNumInnerAnalyticQStrat(7)
-tdacusticsl = @discretise SL[j′,j] == -1.0E[j′] j∈V j′∈W
+tdacusticsl = @discretise SL[j′,j] == -1.0e[j′] j∈V j′∈W
xacusticsl = solve(tdacusticsl)
-#corregere da qui
-import Plots
-Plots.plot(xefie[1,:])
-
-import Plotly
-fcr, geo = facecurrents(xefie[:,125], X)
-Plotly.plot(patch(geo, norm.(fcr)))
-
-
-
-
+Z = assemble(SL,W,V)
+#corregere da qui
+#import Plots
+#Plots.plot(xacusticsl[1,250:300])
+
+#import Plotly
+#fcr, geo = facecurrents(xefie[:,125], X)
+#Plotly.plot(patch(geo, norm.(fcr)))
+
+import BEAST.ConvolutionOperators
+
+for a in 0:300
+ za=ConvolutionOperators.timeslice(Z,a)
+ #zawilton=ConvolutionOperators.timeslice(Zs,a)
+ for i in 1:numfunctions(X)
+ for j in 1:numfunctions(X)
+ if (za[i,j])<0.0
+ println(a,[i,j]," ",za[i,j]," ",)
+ end
+ end
+ end
+end
+
+name=readline()
+parse(Float64, name)
Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
ω = collect(ω0 .+ (0:Nt-1)*Δω)
@@ -48,4 +59,5 @@ _, i1 = findmin(abs.(ω.-1.0))
ue = Xefie[:,i1] / fouriertransform(gaussian)(ω1)
fcr, geo = facecurrents(ue, X)
-Plotly.plot(patch(geo, norm.(fcr)))
\ No newline at end of file
+Plotly.plot(patch(geo, norm.(fcr)))
+
diff --git a/src/helmholtz3d/timedomain/tdhh3dexc.jl b/src/helmholtz3d/timedomain/tdhh3dexc.jl
index 045ea0b8..9119ca5c 100644
--- a/src/helmholtz3d/timedomain/tdhh3dexc.jl
+++ b/src/helmholtz3d/timedomain/tdhh3dexc.jl
@@ -10,6 +10,7 @@ function planewave(direction, speedoflight::Number, signature, amplitude=one(spe
PlaneWaveHH3DTD(direction, speedoflight, signature, amplitude)
end
+scalartype(p::PlaneWaveHH3DTD) = scalartype(p.amplitude)
*(a, f::PlaneWaveHH3DTD) = PlaneWaveHH3DTD(f.direction, f.speed_of_light, f.signature, a*f.amplitude)
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
index 9bbbcca0..d23edece 100644
--- a/src/maxwell/timedomain/acustictdops.jl
+++ b/src/maxwell/timedomain/acustictdops.jl
@@ -31,8 +31,8 @@ end #of the module
export TDAcustic3D
-defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1)
-
+defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = #=nothing=# AllAnalyticalQStrat(1)
+#nothing goes in hybrid qr, allanalytical goes in zuccottirule
@@ -44,24 +44,96 @@ end
quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
p, testel, q, trialel, r, timeel, qd, ::AllAnalyticalQStrat) = ZuccottiRule(1.0)
+
-function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- testels, trialels, timeels, quadstrat::OuterNumInnerAnalyticQStrat)
+
+function quaddata(operator::AcusticSingleLayerTDIO,
+ test_local_space, trial_local_space, time_local_space,
+ test_element, trial_element, time_element, quadstrat::Nothing)
- dmax = numfunctions(timerefs)-1
- bn = binomial.((0:dmax),(0:dmax)')
+ dmax = numfunctions(time_local_space)-1
+ bn = binomial.((0:dmax),(0:dmax)')
- V = eltype(testels[1].vertices)
- ws = WiltonInts84.workspace(V)
- # quadpoints(testrefs, testels, (3,)), bn, ws
- quadpoints(testrefs, testels, (quadstrat.outer_rule,)), bn, ws
-end
+ V = eltype(test_element[1].vertices)
+ ws = WiltonInts84.workspace(V)
+ order = 4
+ @show order
+ quadpoints(test_local_space, test_element, (order,)), bn, ws
+ end
-quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- p, testel, q, trialel, r, timeel, qd, ::OuterNumInnerAnalyticQStrat) = WiltonInts84Strat(qd[1][1,p],qd[2],qd[3])
+
+ # See: ?BEAST.quadrule for help
+ function quadrule(operator::AcusticSingleLayerTDIO,
+ test_local_space, trial_local_space, time_local_space,
+ p, test_element, q, trial_element, r, time_element,
+ quad_data, quadstrat::Nothing)
+
+ # WiltonInts84Strat(quad_data[1,p])
+ qd = quad_data
+ HybridZuccottiWiltonStrat(qd[1][1,p],qd[2],qd[3])
+
+ end
+
+
+ function innerintegrals!(zlocal, operator::AcusticSingleLayerTDIO,
+ test_point,
+ test_local_space, trial_local_space, time_local_space,
+ test_element, trial_element, time_element,
+ quad_rule::HybridZuccottiWiltonStrat, quad_weight)
+
+ # error("Here!!!")
+
+ dx = quad_weight
+ x = cartesian(test_point)
+ # n = normal(test_point)
+
+ # a = trial_element[1]
+ # ξ = x - dot(x -a, n) * n
+
+ r = time_element[1]
+ R = time_element[2]
+ @assert r < R
+
+ N = max(degree(time_local_space), 1)
+ ∫G, ∫vG, ∫∇G = WiltonInts84.wiltonints(
+ trial_element[1],
+ trial_element[2],
+ trial_element[3],
+ x, r, R, Val{2}, quad_rule.workspace)
+
+ a = dx / (4*pi)
+ #D = operator.num_diffs
+ D=0
+ @assert D == 0
+ @assert numfunctions(test_local_space) == 1
+ @assert numfunctions(trial_local_space) == 1
+
+ @inline function tmRoR_sl(d, iG)
+ sgn = isodd(d) ? -1 : 1
+ r = sgn * iG[d+2]
+ end
+
+ # bns = quad_rule.binomials
+
+ @assert D == 0
+ for k in 1 : numfunctions(time_local_space)
+ d = k - 1
+ d < D && continue
+ q = reduce(*, d-D+1:d ,init=1)
+ zlocal[1,1,k] += a * q * tmRoR_sl(d-D, ∫G)
+ end # k
+ end
-function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::ZuccottiStrat) where T
+
+
+qpclineline=10^8#quasiparallel case controller
+qpclinetriang=10^10
+qpctriangtriang=10^13
+
+const qpc=[qpclineline, qpclinetriang, qpctriangtriang]
+
+function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::ZuccottiRule) where T
t1=ι[1]
t2=ι[2]
@@ -92,11 +164,13 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
end
end
end
- #print(τ[1]," ",τ[2]," ",τ[3]," ",σ[1]," ",σ[2]," ",σ[3]," ",t1," ",t2," ",)
+
if hits==3
- z[1,1,1]+=TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2)
+ #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nt1=",t1,"\nt2=",t2)
+
+ z[1,1,1]+=(TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2))/(4*π)
elseif hits==2
- #pay attention with double layer and index permutation
+ #pay attention with double layer and index permutation or in whatever case has nx
if mod1(a1index +1,3)==a2index
a3index=mod1(a1index-1,3)
else
@@ -104,19 +178,94 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
end
if mod1(b1index+1,3)==b2index
- b3index=mod1(a1index-1,3)
+ b3index=mod1(b1index-1,3)
else
- b3index=mod1(a1index+1,3)
+ b3index=mod1(b1index+1,3)
end
+ #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
+
- z[1,1,1]+=TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[a3index],t1,t2)
+ z[1,1,1]+=(TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
elseif hits==1
a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
-
- z[1,1,1]+=TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2)
+ #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
+ z[1,1,1]+=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
else
- z[1,1,1]+=inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2)
+ #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
+ z[1,1,1]+=(inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2,qpc))/(4*π)
end
#for the moment sos=1 but I will correct this
+end
+
+
+function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::HybridZuccottiWiltonStrat) where T
+
+ t1=ι[1]
+ t2=ι[2]
+
+
+ @assert t2 > t1
+
+ sos = op.speed_of_light
+
+ a1index,a2index,a3index,b1index,b2index,b3index=0,0,0,0,0,0
+ Tt = eltype(eltype(τ.vertices))
+ hits = 0
+ dtol = 1.0e3 * eps(Tt)
+ dmin2 = floatmax(Tt)
+ for t in 1:3
+ for s in 1:3
+ d2 = LinearAlgebra.norm_sqr(τ[t]-σ[s])
+ dmin2 = min(dmin2, d2)
+ if (d2 < dtol)
+ hits+=1
+ if hits==1
+ a1index =t
+ b1index = s
+ elseif hits==2
+ a2index=t
+ b2index=s
+ end
+ end
+ end
+ end
+
+ if hits==3
+ #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nt1=",t1,"\nt2=",t2)
+
+ z[1,1,1]+=(TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2))/(4*π)
+ elseif hits==2
+ #pay attention with double layer and index permutation or in whatever case has nx
+ if mod1(a1index +1,3)==a2index
+ a3index=mod1(a1index-1,3)
+ else
+ a3index=mod1(a1index+1,3)
+ end
+
+ if mod1(b1index+1,3)==b2index
+ b3index=mod1(b1index-1,3)
+ else
+ b3index=mod1(b1index+1,3)
+ end
+ #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
+
+
+ z[1,1,1]+=(TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
+ elseif hits==1
+ a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
+ b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
+ #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
+ z[1,1,1]+=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
+ else
+ #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
+ XW = qr.outer_quad_points
+ for p in 1 : length(XW)
+ x = XW[p].point
+ w = XW[p].weight
+ innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end
+
+ end
+ #for the moment sos=1 but I will correct this
end
\ No newline at end of file
diff --git a/src/timedomain/tdintegralop.jl b/src/timedomain/tdintegralop.jl
index 5c6b10c2..1e2637e3 100644
--- a/src/timedomain/tdintegralop.jl
+++ b/src/timedomain/tdintegralop.jl
@@ -293,6 +293,12 @@ struct WiltonInts84Strat{T,V,W}
workspace::W
end
+struct HybridZuccottiWiltonStrat{T,V,W}
+ outer_quad_points::T
+ binomials::V
+ workspace::W
+end
+
struct ZuccottiRule{T}
weight::T #Allanalyticalformula
end
diff --git a/test/runtests.jl b/test/runtests.jl
index cf2ca17a..223922c9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -71,7 +71,8 @@ include("test_tdrhs_scaling.jl")
include("test_td_tensoroperator.jl")
include("test_variational.jl")
-
+include("test_ncrossbdm.jl")
+include("test_curl_lagc0d1_lagc0d2.jl")
try
Pkg.installed("BogaertInts10")
From a1f19cfc3a5ef9bf600b94c6f12c241f4fd58448 Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Wed, 17 Jan 2024 00:39:15 +0000
Subject: [PATCH 06/24] CompatHelper: bump compat for SparseMatrixDicts to 0.2,
(keep existing compat)
---
Project.toml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Project.toml b/Project.toml
index 8c427f99..1b5c0c44 100644
--- a/Project.toml
+++ b/Project.toml
@@ -48,7 +48,7 @@ NestedUnitRanges = "0.2"
Requires = "1"
SauterSchwab3D = "0.1"
SauterSchwabQuadrature = "2.3.0"
-SparseMatrixDicts = "< 0.2.8"
+SparseMatrixDicts = "< 0.2.8, 0.2"
SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1, 2"
StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1"
WiltonInts84 = "0.2.5"
From 06589747517bcd74fc78ebd9b64af11e8711e332 Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Sat, 10 Feb 2024 00:35:12 +0000
Subject: [PATCH 07/24] CompatHelper: bump compat for LinearMaps to 3, (keep
existing compat)
---
Project.toml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Project.toml b/Project.toml
index 8c427f99..6690ed8a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,7 +43,7 @@ FastGaussQuadrature = "0.3, 0.4, 0.5, 1"
FillArrays = "0.11, 0.12, 0.13, 1"
IterativeSolvers = "0.9"
LiftedMaps = "0.5.1"
-LinearMaps = "3.7 - 3.9"
+LinearMaps = "3.7 - 3.9, 3"
NestedUnitRanges = "0.2"
Requires = "1"
SauterSchwab3D = "0.1"
From fb0f1a7847195e37f75fe2e98bfc95c771939590 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Mon, 4 Mar 2024 17:02:05 +0100
Subject: [PATCH 08/24] ncrossbdm permutation + minor changes
---
examples/disabled/tdhh3d_neumann.jl | 10 ++--
examples/tdacusticsinglelayer.jl | 34 +++++++-----
src/bases/local/ncrossbdmlocal.jl | 22 ++++++++
src/maxwell/timedomain/acustictdops.jl | 77 +++++++++++++++++++++-----
4 files changed, 110 insertions(+), 33 deletions(-)
diff --git a/examples/disabled/tdhh3d_neumann.jl b/examples/disabled/tdhh3d_neumann.jl
index 5db0c60b..adf753ab 100644
--- a/examples/disabled/tdhh3d_neumann.jl
+++ b/examples/disabled/tdhh3d_neumann.jl
@@ -2,8 +2,8 @@ using CompScienceMeshes
using BEAST
using LinearAlgebra
-G = meshsphere(1.0, 0.30)
-#G = CompScienceMeshes.meshmobius(h=0.1)
+#G = meshsphere(1.0, 0.30)
+G = CompScienceMeshes.meshmobius(h=0.2)
c = 1.0
S = BEAST.HH3DSingleLayerTDBIO(c)
@@ -37,8 +37,8 @@ P = timebasiscxd0(Δt, Nt)
# assemble the right hand side
bd = assemble(n⋅h, X ⊗ P)
-Z1d = assemble(Id ⊗ Id, X ⊗ P, X ⊗ P, Val{:bandedstorage})
-Z0d = assemble(D, X ⊗ P, X ⊗ P, Val{:bandedstorage})
+Z1d = assemble(Id ⊗ Id, X ⊗ P, X ⊗ P)
+Z0d = assemble(D, X ⊗ P, X ⊗ P)
Zd = Z0d + (-0.5)*Z1d
u = marchonintime(inv(Zd[:,:,1]), Zd, bd, Nt)
@@ -47,7 +47,7 @@ Zs = assemble(S, X ⊗ δ, X ⊗ P)
v = marchonintime(inv(Zs[:,:,1]), Zs, -bs, Nt)
-tdacusticsl = @discretise S[j′,j] == -1.0e[j′] j∈ (X ⊗ P) j′∈ (X ⊗ δ)
+tdacusticsl = @discretise D[j′,j] == 1.0(n⋅h)[j′] j∈ (X ⊗ P) j′∈ (X ⊗ δ)
xacusticsl = solve(tdacusticsl)
diff --git a/examples/tdacusticsinglelayer.jl b/examples/tdacusticsinglelayer.jl
index 4f049814..b952de14 100644
--- a/examples/tdacusticsinglelayer.jl
+++ b/examples/tdacusticsinglelayer.jl
@@ -1,11 +1,10 @@
using CompScienceMeshes, BEAST, LinearAlgebra
#Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
-Γ = meshsphere(radius=1.0, h=0.30)
+Γ = CompScienceMeshes.meshsphere(1.0,0.3)#CompScienceMeshes.meshcuboid(1.0,1.0,1.0,0.3)
X = lagrangecxd0(Γ)
-Δt = 0.16
-Nt = 300
+Δt, Nt = 0.08, 600
T = timebasisshiftedlagrange(Δt, Nt, 0)
U = timebasisdelta(Δt, Nt)
@@ -19,7 +18,7 @@ e = BEAST.planewave(point(0,0,1), 1.0, gaussian)
@hilbertspace j
@hilbertspace j′
-SL = TDAcustic3D.acusticsinglelayer(speedofsound=1.0, numdiffs=1)
+SL = TDAcustic3D.acusticsinglelayer(speedofsound=1.0, numdiffs=0)
# BEAST.@defaultquadstrat (SL, W, V) BEAST.OuterNumInnerAnalyticQStrat(7)
tdacusticsl = @discretise SL[j′,j] == -1.0e[j′] j∈V j′∈W
@@ -27,27 +26,34 @@ xacusticsl = solve(tdacusticsl)
Z = assemble(SL,W,V)
#corregere da qui
-#import Plots
-#Plots.plot(xacusticsl[1,250:300])
+import Plots
-#import Plotly
+Plots.plot(xacusticsl[1,1:300],label="Current")
+Plots.xlabel!("t")
+Plots.savefig("stablecurrent.png")
+
+xacusticsl[1,200:300]
+#pval=ConvolutionOperators.polyvals(Z)
+import Plotly
#fcr, geo = facecurrents(xefie[:,125], X)
#Plotly.plot(patch(geo, norm.(fcr)))
import BEAST.ConvolutionOperators
-for a in 0:300
+for a in 1:1
za=ConvolutionOperators.timeslice(Z,a)
- #zawilton=ConvolutionOperators.timeslice(Zs,a)
- for i in 1:numfunctions(X)
- for j in 1:numfunctions(X)
- if (za[i,j])<0.0
- println(a,[i,j]," ",za[i,j]," ",)
+ zawilton=ConvolutionOperators.timeslice(Zs,a)
+ for i in 1:10 #numfunctions(X)
+ #for j in 1:numfunctions(X)
+ if norm(za[i,i]-zawilton[i,i])>=10^(-6)
+ println(a," ",i," ",za[i,i]," ",zawilton[i,i])
end
- end
+ #end
end
end
+
+
name=readline()
parse(Float64, name)
diff --git a/src/bases/local/ncrossbdmlocal.jl b/src/bases/local/ncrossbdmlocal.jl
index b7dc7504..fe2a2a85 100644
--- a/src/bases/local/ncrossbdmlocal.jl
+++ b/src/bases/local/ncrossbdmlocal.jl
@@ -18,3 +18,25 @@ function (f::NCrossBDMRefSpace{T})(p) where T
(value= n × (u*tu)/j, curl=d),
(value= n × (v*tv)/j, curl=d),]
end
+
+const _vert_perms_ncrossbdm = [
+ (1,2,3),
+ (2,3,1),
+ (3,1,2),
+ (2,1,3),
+ (1,3,2),
+ (3,2,1),
+]
+const _dof_perms_ncrossbdm = [
+ (1,2,3,4,5,6),
+ (5,6,1,2,3,4),
+ (3,4,5,6,1,2),
+ (4,3,2,1,6,5),
+ (2,1,6,5,4,3),
+ (6,5,4,3,2,1),
+]
+
+function dof_permutation(::NCrossBDMRefSpace, vert_permutation)
+ i = findfirst(==(tuple(vert_permutation...)), _vert_perms_ncrossbdm)
+ return _dof_perms_ncrossbdm[i]
+end
\ No newline at end of file
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
index d23edece..4a865d87 100644
--- a/src/maxwell/timedomain/acustictdops.jl
+++ b/src/maxwell/timedomain/acustictdops.jl
@@ -31,7 +31,7 @@ end #of the module
export TDAcustic3D
-defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = #=nothing=# AllAnalyticalQStrat(1)
+defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = nothing # AllAnalyticalQStrat(1)
#nothing goes in hybrid qr, allanalytical goes in zuccottirule
@@ -233,10 +233,26 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
if hits==3
#print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nt1=",t1,"\nt2=",t2)
-
- z[1,1,1]+=(TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2))/(4*π)
+ entry=(TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2))/(4*π)
+ if norm(entry)>10.0 || entry<-0.0001
+ println("quadrature ",τ[1]," ",τ[2]," ",τ[3]," ", t1," ",t2," ",entry)
+ XW = qr.outer_quad_points
+ for p in 1 : length(XW)
+ x = XW[p].point
+ w = XW[p].weight
+ innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end
+ else
+ z[1,1,1]+=max(entry,0.0)
+ end
elseif hits==2
- #pay attention with double layer and index permutation or in whatever case has nx
+ #= XW = qr.outer_quad_points
+ for p in 1 : length(XW)
+ x = XW[p].point
+ w = XW[p].weight
+ innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end=#
+ #pay attention with double layer and index permutation or whatever case has n cross product
if mod1(a1index +1,3)==a2index
a3index=mod1(a1index-1,3)
else
@@ -250,22 +266,55 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
end
#print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
-
- z[1,1,1]+=(TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
+ entry=(TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
+ if norm(entry)>10.0 || entry<-0.0001
+ println("quadrature ",τ[a1index]," ",τ[a2index]," ",τ[a3index]," ",σ[b1index]," ",σ[b2index]," ",σ[b3index], t1," ",t2," ",entry)
+ XW = qr.outer_quad_points
+ for p in 1 : length(XW)
+ x = XW[p].point
+ w = XW[p].weight
+ innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end
+ else
+ z[1,1,1]+=max(0.0,entry)
+ end
elseif hits==1
- a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
- b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
- #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
- z[1,1,1]+=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
- else
- #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
- XW = qr.outer_quad_points
+ #=XW = qr.outer_quad_points
for p in 1 : length(XW)
x = XW[p].point
w = XW[p].weight
innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end=#
+ a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
+ b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
+ #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
+ entry=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
+ if norm(entry)>10.0 || entry<-0.0001
+ println("quadrature ",τ[a1index]," ",τ[a2index]," ",τ[a3index]," ",σ[b2index]," ",σ[b3index], t1," ",t2," ",entry)
+
+ XW = qr.outer_quad_points
+ for p in 1 : length(XW)
+ x = XW[p].point
+ w = XW[p].weight
+ innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end
+ else
+ z[1,1,1]+=max(0.0,entry)
+ end
+ else
+ entry=0.1#(inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2,qpc))/(4*π)
+ #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
+ if norm(entry)>0.0 || entry<-0.0001
+ XW = qr.outer_quad_points
+ for p in 1 : length(XW)
+ x = XW[p].point
+ w = XW[p].weight
+ innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
+ end
+
+ else
+ z[1,1,1]+=max(0.0,entry)
end
-
end
#for the moment sos=1 but I will correct this
end
\ No newline at end of file
From 7d6707e7743af615f24a92307c65e6fb1fbb9245 Mon Sep 17 00:00:00 2001
From: CompatHelper Julia <compathelper_noreply@julialang.org>
Date: Fri, 10 May 2024 00:38:56 +0000
Subject: [PATCH 09/24] CompatHelper: bump compat for CompScienceMeshes to 0.7,
(keep existing compat)
---
Project.toml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Project.toml b/Project.toml
index 6f052b4e..3108e7eb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -35,7 +35,7 @@ AbstractTrees = "0.4.4"
BlockArrays = "0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
CollisionDetection = "0.1.5"
Combinatorics = "0.7, 1"
-CompScienceMeshes = "0.6.0"
+CompScienceMeshes = "0.6.0, 0.7"
Compat = "2, 3, 4"
ConvolutionOperators = "0.4"
ExtendableSparse = "1.4"
From a75e2d0cec6d632edb128008484f04176b535576 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 15 May 2024 14:41:14 +0200
Subject: [PATCH 10/24] minor modifications
---
examples/tdacusticsinglelayer.jl | 7 ++++---
src/maxwell/timedomain/acustictdops.jl | 7 ++++---
2 files changed, 8 insertions(+), 6 deletions(-)
diff --git a/examples/tdacusticsinglelayer.jl b/examples/tdacusticsinglelayer.jl
index b952de14..a55aa005 100644
--- a/examples/tdacusticsinglelayer.jl
+++ b/examples/tdacusticsinglelayer.jl
@@ -1,10 +1,10 @@
using CompScienceMeshes, BEAST, LinearAlgebra
#Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
-Γ = CompScienceMeshes.meshsphere(1.0,0.3)#CompScienceMeshes.meshcuboid(1.0,1.0,1.0,0.3)
+Γ = CompScienceMeshes.meshsphere(1.0,0.4)#CompScienceMeshes.meshcuboid(1.0,1.0,1.0,0.3)
X = lagrangecxd0(Γ)
-Δt, Nt = 0.08, 600
+Δt, Nt = 0.5, 200
T = timebasisshiftedlagrange(Δt, Nt, 0)
U = timebasisdelta(Δt, Nt)
@@ -25,10 +25,11 @@ tdacusticsl = @discretise SL[j′,j] == -1.0e[j′] j∈V j′∈W
xacusticsl = solve(tdacusticsl)
Z = assemble(SL,W,V)
+
#corregere da qui
import Plots
-Plots.plot(xacusticsl[1,1:300],label="Current")
+Plots.plot(xacusticsl[1,1:300],label="Current",ylim=[-0.0001,0.0001])
Plots.xlabel!("t")
Plots.savefig("stablecurrent.png")
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
index 4a865d87..f82e12bd 100644
--- a/src/maxwell/timedomain/acustictdops.jl
+++ b/src/maxwell/timedomain/acustictdops.jl
@@ -31,7 +31,7 @@ end #of the module
export TDAcustic3D
-defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = nothing # AllAnalyticalQStrat(1)
+defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1)
#nothing goes in hybrid qr, allanalytical goes in zuccottirule
@@ -189,8 +189,9 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
elseif hits==1
a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
- #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
+ print("\na1=",τ[a1index],"\na2=",τ[a2index],"\na3=",τ[a3index],"\nb1=",σ[b1index],"\nb2=",σ[b2index],"\nb3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
z[1,1,1]+=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
+
else
#print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
z[1,1,1]+=(inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2,qpc))/(4*π)
@@ -290,7 +291,7 @@ function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}
#print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
entry=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
if norm(entry)>10.0 || entry<-0.0001
- println("quadrature ",τ[a1index]," ",τ[a2index]," ",τ[a3index]," ",σ[b2index]," ",σ[b3index], t1," ",t2," ",entry)
+ println("quadrature ",τ[a1index]," ",τ[a2index]," ",τ[a3index]," ",σ[b2index]," ",σ[b3index], t1," ",t2," ",entry, " ",)
XW = qr.outer_quad_points
for p in 1 : length(XW)
From d204aac8efc53996252566eed01358f04eca4509 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Thu, 13 Jun 2024 16:20:40 +0200
Subject: [PATCH 11/24] Building new structure for all analytical ints
---
Project.toml | 1 +
src/maxwell/timedomain/acustictdops.jl | 6 +-
src/timedomain/analyticalints.jl | 218 +++++++++++++++++++++++++
3 files changed, 221 insertions(+), 4 deletions(-)
create mode 100644 src/timedomain/analyticalints.jl
diff --git a/Project.toml b/Project.toml
index 0f162a22..cc4444bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -26,6 +26,7 @@ SauterSchwab3D = "0a13313b-1c00-422e-8263-562364ed9544"
SauterSchwabQuadrature = "535c7bfe-2023-5c1d-b712-654ef9d93a38"
SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseMatrixDicts = "5cb6c4b0-9b79-11e8-24c9-f9621d252589"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
TimeDomainBEMInt = "0303f91b-5811-4998-82cd-e0687f55e8f9"
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
index f82e12bd..e25c5078 100644
--- a/src/maxwell/timedomain/acustictdops.jl
+++ b/src/maxwell/timedomain/acustictdops.jl
@@ -31,14 +31,12 @@ end #of the module
export TDAcustic3D
-defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1)
+defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1) #nothing
#nothing goes in hybrid qr, allanalytical goes in zuccottirule
-
-
function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
+ testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
return nothing
end
diff --git a/src/timedomain/analyticalints.jl b/src/timedomain/analyticalints.jl
new file mode 100644
index 00000000..bce9909e
--- /dev/null
+++ b/src/timedomain/analyticalints.jl
@@ -0,0 +1,218 @@
+function minmax1d(vertex,edge)
+ T = eltype(τ[1])
+ m = norm(vertex-edge[1])
+ M = m
+ s=edge[1]-edge[2]
+ s/=norm(s)
+ ev1=edge[1]-vertex
+ x0=(edge[1]-dot(ev1,s)*s)
+ a=(edge[2]-x0)*s
+ b=(edge[1]-x0)*s
+ if a<=0 && b>=0
+ m=norm(vertex-x0)
+ abs(a)<abs(b) && (M=norm(vertex-edge[2]))
+
+ else
+ for j in 1:length(edge)
+ q = edge[j]
+ d = norm(vertex-q)
+ d < m && (m=d)
+ d > M && (M=d)
+ end
+ end
+ return m, M
+end
+
+function rings1d(τ, σ, ΔR)
+ m, M = minmaxdist(τ, σ)
+ r0 = floor(Int, m/ΔR) + 1
+ r1 = ceil(Int, M/ΔR+1)
+ r0 : r1
+end
+
+
+
+function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
+ testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
+
+ dimU=dimension(testels)
+ dimV=dimension(trialels)
+ #rigerenerare delta R
+
+ if dimU+dimV==1
+ #testelsboundary=skeleton(testels,dimU-1)
+ # trialelsboundary=skeleton(trialels,dimV-1)
+
+ numnodes=length(testels)
+ numedges=length(trialels)
+
+ datavertexedge=Array{TimeDomainBEMInt.edgevertexgeo{T,P}, 2}(undef, numnodes, numedges)
+ rings=Array{UnitRange{Int},2}(undef, numnodes, numedges)
+ datarings=Array{Vector{Tuple{Int,Vector}},2}(undef,numnodes,numedges)#il type va bene
+
+ #fill datarings with zeross !
+
+ for p in 1:numnodes
+ τ = chart(testels,p)#testels[p]
+ for q in 1:numedges
+ σ = chart(trialels,q)
+ edgevertgeo=TimeDomainBEMInt.edgevertexinteraction(τ,σ[1],σ[2])
+ datavertexedge[p,q]=edgevertgeo
+ a,b=edgevertegeo.extint0[1],edgevertegeo.extint0[2]
+ rngs=rings1d(τ,σ,ΔR)
+ rings[p,q]=rngs
+ datarings[p,q]=[0,[0.0,0.0]]
+ for r in rngs
+ r > numfunctions(timebasisfunction) && continue #serve?
+ ι = ring(r,ΔR)
+
+ # compute interactions between reference shape functions
+ #fill!(z, 0)
+ rp=τ #se e un simplex ok se no va messo chart(τ,1).vertices credo
+ t2=ι[2]#needs a check
+ extint=TimeDomainBEMInt.edgevertexinteraction(t2,edgevertgeo)
+ push!(datarings[p,q],extint)
+
+ # qr = quadrule(op, U, V, W, p, τ, q, σ, r, ι, qd, quadstrat)
+ #momintegrals!(z, op, U, V, W, τ, σ, ι, qr)
+ end
+ end
+ end
+
+ return datavertexedge
+ else
+ return "devo ancora scrivere"
+ end
+end
+
+
+nnodes=length(nodes)
+edge1=chart(edges,p)
+edge2=chart(edges,q)
+cnnct=connectivity(edges,nodes)
+vertind1=cnnct[1:nnodes,p].nzind
+vertsgn1=cnnct[1:nnodes,p].nzval
+vertind2=cnnct[1:nnodes,q].nzind
+vertsgn2=cnnct[1:nnodes,q].nzval
+
+if vertsgn1[1]==1
+ a1,a2=edge1[1],edge1[2]
+else
+ a2,a1=edge1[1],edge1[2]
+end
+
+if vertsgn2[1]==1
+ b1,b2=edge1[1],edge1[2]
+else
+ b2,b1=edge1[1],edge1[2]
+end
+
+geo1,rings1,datarings1=edgevertexgeo[vertind1[1],q],rings[vertind1[1],q],datarings[vertind1[1],q]
+geo2,rings2,datarings2=edgevertexgeo[vertind1[2],q],rings[vertind1[2],q],datarings[vertind1[2],q]
+geo3,rings3,datarings3=edgevertexgeo[vertind2[1],p],rings[vertind2[1],p],datarings[vertind2[1],p]
+geo4,rings4,datarings4=edgevertexgeo[vertind2[2],p],rings[vertind2[2],p],datarings[vertind2[2],p]
+
+geo=[geo1,geo2,geo3,geo4]
+rings=[rings1,rings2,rings3,rings4]
+datarings=[datarings1,datarings2,datarings3,datarings4]
+
+
+
+function intlinelineglobal(a1,a2,b1,b2,geo,rings,datatimes,parcontrol,UB::Type{Val{N}}) where N
+
+ #nedges=length(edges)
+
+
+ #vertices=[a1,a2,a1′,a2′]
+ vertices=[a1,a2,b1,b2]
+ l12=norm(vertices[1]-vertices[2])
+ l12′=norm(vertices[3]-vertices[4])
+
+
+ #geo1=edgevertexinteraction(a1,a1′,a2′)
+ #geo2=edgevertexinteraction(a2,a1′,a2′)
+ #geo3=edgevertexinteraction(a1′,a1,a2)
+ #geo4=edgevertexinteraction(a2′,a1,a2)
+
+ #datatime=Array{Tuple}(undef,2,4)?
+ I = maketuple(eltype(a1), UB)
+ K = maketuple(typeof(a1), UB)
+
+ x=geo[3].tangent
+
+ #z=cross(a12′,x)
+ #J=norm(z)
+ #if J blablabla
+ #z /= J
+ #h=dot(r22,z)
+ hdir=cross(geo[1].tangent,x)
+ n=hdir/norm(hdir)
+ sgnn=[+1,-1,-1,+1]
+ h=dot(a2-b2,n)
+ sgnh=[+1,-1,+1,-1]
+ angletot=0.0
+ dminv=Vector{eltype(edge1[1])}(undef, 4)
+ ξ=Vector{typeof(edge1[1])}(undef, 4)
+ for j in 1:4
+ dminv[j]=geo[j].dmin
+ dmaxv[j]=geo[j].dmax
+ v=vertices[j]
+ ξ[j]=v-n*h*sgnh[j]*sgnn[j]
+ angletot+=anglecontribution(ξ[j],sgnn[j]*n,geo[j])
+ end
+ if abs(angletot-2π)<100*eps(eltype(edge1[1]))
+ dmin=abs(h)
+ else
+ dmin=min(dminv[1],dminv[2],dminv[3],dminv[4])
+ end
+
+ dmax=max(dmaxv[1],dmaxv[2],dmaxv[3],dmaxv[4])
+
+ r0 = floor(Int, dmin/ΔR) + 1
+ r1 = ceil(Int, dmax/ΔR+1) #recuperare deltaR
+ ringtot = r0 : r1
+
+ allint=Vector{typeof((I,K))}(undef,r1-r0+2)
+ fill!(allint,(I,K))
+ if norm(hdir) < (parcontrol[1])*eps(typeof(temp1))
+ I=intparallelsegment(a1,a2,b1,b2,temp1,temp2)[1] #attenzione qui non compatibile con quello che stiamo scrivendo
+ else
+ n=hdir/norm(hdir)
+ sgnn=[+1,-1,-1,+1]
+ h=dot(a2-a2′,n)
+ sgnh=[+1,-1,+1,-1]
+ for j in 1:4
+ for i in ringtot[1]:(relrings[j][1]-1)
+
+ P,Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],0,[0,0],i*ΔR,UB)
+ allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
+ allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
+ P,Q = arcsegcontribution(v,ξ[j],-sgnn[j]*n,sgnh[j]*h,geo[j],0,[0,0],(i-1)*ΔR,UB)
+ allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
+ allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
+ end
+ for i in relrings[j]
+
+ #shall I put some check like i*deltaR > h
+ P, Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-relrings[j][1]+2][1],datarings[j][i-relrings[j][1]+2][2],i*ΔR,UB) #salviamo
+ allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
+ allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
+ P, Q = arcsegcontribution(v,ξ[j],-sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-relrings[j][1]+1][1],datarings[j][i-relrings[j][1]+1][2],(i-1)*ΔR,UB)
+ allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
+ allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
+ end #probabilmente va bene cosi anche con ceil(int,frac) invece di ceil(int,frac+1)
+ allint[i-ringtot[1]+2][1]=multiply(allint[i-ringtot[1]+2][1],1/(l12′*l12*norm(hdir)))
+ allint[i-ringtot[1]+2][2]=multiply(allint[i-ringtot[1]+2][2],1/(l12′*l12*norm(hdir)))
+ #= for i in (relrings[j][2]+1):ringtot[j][2]
+
+ #shall I put some check like i*deltaR > h
+ saveP[i],saveQ[i] = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],1,[a,b],i*ΔR,UB) #dadefinire a e b
+ save,and,subtract = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],1,[a,b],(i-1)*ΔR,UB)
+ end =# #sembra che non serva a causa di ceil(int,frac+1)!
+ end
+ #I+=(1/abs(r12′[2]*r12[1]))*(3*(temp2^2-temp1^2)*d[1]-2*(temp2^3-temp1^3)*d[2])
+ end
+
+ return ringtot,allint #missing buidgrad since it is not yet adapted for int line line
+end
+
From c2d4539432b011e8054daa8d0d608d19e704749a Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Mon, 17 Jun 2024 17:10:03 +0200
Subject: [PATCH 12/24] Update analyticalints.jl
---
src/timedomain/analyticalints.jl | 88 ++++++++++++++++++--------------
1 file changed, 51 insertions(+), 37 deletions(-)
diff --git a/src/timedomain/analyticalints.jl b/src/timedomain/analyticalints.jl
index bce9909e..97e4359a 100644
--- a/src/timedomain/analyticalints.jl
+++ b/src/timedomain/analyticalints.jl
@@ -33,13 +33,13 @@ end
function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
+ testels::Vector{Simplex{3,0,3,1,T}}, trialels::Vector{Simplex{3,1,1,2,T}}, timeels, quadstrat::AllAnalyticalQStrat) where T
dimU=dimension(testels)
dimV=dimension(trialels)
#rigerenerare delta R
-
- if dimU+dimV==1
+#quaddata 1D
+ #if dimU+dimV==1
#testelsboundary=skeleton(testels,dimU-1)
# trialelsboundary=skeleton(trialels,dimV-1)
@@ -80,45 +80,59 @@ function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
end
return datavertexedge
- else
- return "devo ancora scrivere"
- end
+ #else
+ # return "devo ancora scrivere"
+ #end
end
+function quaddata2Dee(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
+ testels::Vector{Simplex{3,1,1,2,T}}, trialels::Vector{Simplex{3,1,1,2,T}}, timeels, quadstrat::AllAnalyticalQStrat)
+ nnodes=length(nodes)
-nnodes=length(nodes)
-edge1=chart(edges,p)
-edge2=chart(edges,q)
-cnnct=connectivity(edges,nodes)
-vertind1=cnnct[1:nnodes,p].nzind
-vertsgn1=cnnct[1:nnodes,p].nzval
-vertind2=cnnct[1:nnodes,q].nzind
-vertsgn2=cnnct[1:nnodes,q].nzval
-
-if vertsgn1[1]==1
- a1,a2=edge1[1],edge1[2]
-else
- a2,a1=edge1[1],edge1[2]
-end
+
+ totrings=Array{UnitRange{Int},2}(undef, numedges, numedges)
+ datavalues=Array{Vector{Tuple{Int,Vector}},2}(undef,numedges,numedges)
+ cnnct=connectivity(edges,nodes)
+
+ for p in 1:numdges
+ for q in 1:numedges
+ edge1=chart(edges,p)
+ edge2=chart(edges,q)
+
+ vertind1=cnnct[1:nnodes,p].nzind
+ vertsgn1=cnnct[1:nnodes,p].nzval
+ vertind2=cnnct[1:nnodes,q].nzind
+ vertsgn2=cnnct[1:nnodes,q].nzval
+
+ if vertsgn1[1]==1
+ a1,a2=edge1[1],edge1[2]
+ else
+ a2,a1=edge1[1],edge1[2]
+ end
-if vertsgn2[1]==1
- b1,b2=edge1[1],edge1[2]
-else
- b2,b1=edge1[1],edge1[2]
-end
+ if vertsgn2[1]==1
+ b1,b2=edge2[1],edge2[2]
+ else
+ b2,b1=edge2[1],edge2[2]
+ end
-geo1,rings1,datarings1=edgevertexgeo[vertind1[1],q],rings[vertind1[1],q],datarings[vertind1[1],q]
-geo2,rings2,datarings2=edgevertexgeo[vertind1[2],q],rings[vertind1[2],q],datarings[vertind1[2],q]
-geo3,rings3,datarings3=edgevertexgeo[vertind2[1],p],rings[vertind2[1],p],datarings[vertind2[1],p]
-geo4,rings4,datarings4=edgevertexgeo[vertind2[2],p],rings[vertind2[2],p],datarings[vertind2[2],p]
+ geo1,rings1,datarings1=edgevertexgeo[vertind1[1],q],rings[vertind1[1],q],datarings[vertind1[1],q]
+ geo2,rings2,datarings2=edgevertexgeo[vertind1[2],q],rings[vertind1[2],q],datarings[vertind1[2],q]
+ geo3,rings3,datarings3=edgevertexgeo[vertind2[1],p],rings[vertind2[1],p],datarings[vertind2[1],p]
+ geo4,rings4,datarings4=edgevertexgeo[vertind2[2],p],rings[vertind2[2],p],datarings[vertind2[2],p]
-geo=[geo1,geo2,geo3,geo4]
-rings=[rings1,rings2,rings3,rings4]
-datarings=[datarings1,datarings2,datarings3,datarings4]
+ geo=[geo1,geo2,geo3,geo4]
+ rings=[rings1,rings2,rings3,rings4]
+ datarings=[datarings1,datarings2,datarings3,datarings4]
+ totrings[p,q],datavalues[p,q]=intlinelineglobal(a1,a2,b1,b2,geo,rings,datarings,[10^6,10^6,10^6],Val{0})
+ end
+ end
+ return totrings,datavalues
+end
-function intlinelineglobal(a1,a2,b1,b2,geo,rings,datatimes,parcontrol,UB::Type{Val{N}}) where N
+function intlinelineglobal(a1,a2,b1,b2,geo,rings,datarings,parcontrol,UB::Type{Val{N}}) where N
#nedges=length(edges)
@@ -182,7 +196,7 @@ function intlinelineglobal(a1,a2,b1,b2,geo,rings,datatimes,parcontrol,UB::Type{V
h=dot(a2-a2′,n)
sgnh=[+1,-1,+1,-1]
for j in 1:4
- for i in ringtot[1]:(relrings[j][1]-1)
+ for i in ringtot[1]:(rings[j][1]-1)
P,Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],0,[0,0],i*ΔR,UB)
allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
@@ -191,13 +205,13 @@ function intlinelineglobal(a1,a2,b1,b2,geo,rings,datatimes,parcontrol,UB::Type{V
allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
end
- for i in relrings[j]
+ for i in rings[j]
#shall I put some check like i*deltaR > h
- P, Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-relrings[j][1]+2][1],datarings[j][i-relrings[j][1]+2][2],i*ΔR,UB) #salviamo
+ P, Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-rings[j][1]+2][1],datarings[j][i-rings[j][1]+2][2],i*ΔR,UB) #salviamo
allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
- P, Q = arcsegcontribution(v,ξ[j],-sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-relrings[j][1]+1][1],datarings[j][i-relrings[j][1]+1][2],(i-1)*ΔR,UB)
+ P, Q = arcsegcontribution(v,ξ[j],-sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-rings[j][1]+1][1],datarings[j][i-rings[j][1]+1][2],(i-1)*ΔR,UB)
allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
end #probabilmente va bene cosi anche con ceil(int,frac) invece di ceil(int,frac+1)
From 68804cb5c5df02f78c27ccf933b70a40f9492177 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 9 Oct 2024 12:47:43 +0200
Subject: [PATCH 13/24] Update runtests.jl
---
test/runtests.jl | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/test/runtests.jl b/test/runtests.jl
index ac4f15cb..b642a5d2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -78,7 +78,8 @@ include("test_variational.jl")
include("test_handlers.jl")
include("test_ncrossbdm.jl")
-include("test_curl_lagc0d1_lagc0d2.jl")include("test_gridfunction.jl")
+include("test_curl_lagc0d1_lagc0d2.jl")
+include("test_gridfunction.jl")
using TestItemRunner
@run_package_tests
From 97f4c1d270894bb215537a4840e6a3b49c84f0c1 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 9 Oct 2024 13:18:40 +0200
Subject: [PATCH 14/24] Cleaning master branch
I moved my changes for time domain integration in a new master called time_domain+integration
---
Project.toml | 2 -
examples/tdacusticsinglelayer.jl | 70 ------
src/BEAST.jl | 1 -
src/maxwell/timedomain/acustictdops.jl | 319 -------------------------
src/quadrature/quadstrats.jl | 4 +-
src/timedomain/analyticalints.jl | 232 ------------------
src/timedomain/tdintegralop.jl | 10 -
7 files changed, 1 insertion(+), 637 deletions(-)
delete mode 100644 examples/tdacusticsinglelayer.jl
delete mode 100644 src/maxwell/timedomain/acustictdops.jl
delete mode 100644 src/timedomain/analyticalints.jl
diff --git a/Project.toml b/Project.toml
index 6b722b59..98f53e4e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,10 +27,8 @@ SauterSchwab3D = "0a13313b-1c00-422e-8263-562364ed9544"
SauterSchwabQuadrature = "535c7bfe-2023-5c1d-b712-654ef9d93a38"
SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseMatrixDicts = "5cb6c4b0-9b79-11e8-24c9-f9621d252589"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-TimeDomainBEMInt = "0303f91b-5811-4998-82cd-e0687f55e8f9"
TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
WiltonInts84 = "a3e2863e-c0ee-5ff6-a523-307a4cdc8724"
diff --git a/examples/tdacusticsinglelayer.jl b/examples/tdacusticsinglelayer.jl
deleted file mode 100644
index a55aa005..00000000
--- a/examples/tdacusticsinglelayer.jl
+++ /dev/null
@@ -1,70 +0,0 @@
-using CompScienceMeshes, BEAST, LinearAlgebra
-#Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
-Γ = CompScienceMeshes.meshsphere(1.0,0.4)#CompScienceMeshes.meshcuboid(1.0,1.0,1.0,0.3)
-
-X = lagrangecxd0(Γ)
-
-Δt, Nt = 0.5, 200
-T = timebasisshiftedlagrange(Δt, Nt, 0)
-U = timebasisdelta(Δt, Nt)
-
-V = X ⊗ T
-W = X ⊗ U
-
-width, delay, scaling = 16.0, 24.0, 1.0
-gaussian = creategaussian(width, delay, scaling)
-e = BEAST.planewave(point(0,0,1), 1.0, gaussian)
-
-@hilbertspace j
-@hilbertspace j′
-
-SL = TDAcustic3D.acusticsinglelayer(speedofsound=1.0, numdiffs=0)
-# BEAST.@defaultquadstrat (SL, W, V) BEAST.OuterNumInnerAnalyticQStrat(7)
-
-tdacusticsl = @discretise SL[j′,j] == -1.0e[j′] j∈V j′∈W
-xacusticsl = solve(tdacusticsl)
-
-Z = assemble(SL,W,V)
-
-#corregere da qui
-import Plots
-
-Plots.plot(xacusticsl[1,1:300],label="Current",ylim=[-0.0001,0.0001])
-Plots.xlabel!("t")
-Plots.savefig("stablecurrent.png")
-
-xacusticsl[1,200:300]
-#pval=ConvolutionOperators.polyvals(Z)
-import Plotly
-#fcr, geo = facecurrents(xefie[:,125], X)
-#Plotly.plot(patch(geo, norm.(fcr)))
-
-import BEAST.ConvolutionOperators
-
-for a in 1:1
- za=ConvolutionOperators.timeslice(Z,a)
- zawilton=ConvolutionOperators.timeslice(Zs,a)
- for i in 1:10 #numfunctions(X)
- #for j in 1:numfunctions(X)
- if norm(za[i,i]-zawilton[i,i])>=10^(-6)
- println(a," ",i," ",za[i,i]," ",zawilton[i,i])
- end
- #end
- end
-end
-
-
-
-name=readline()
-parse(Float64, name)
-
-Xefie, Δω, ω0 = fouriertransform(xefie, Δt, 0.0, 2)
-ω = collect(ω0 .+ (0:Nt-1)*Δω)
-_, i1 = findmin(abs.(ω.-1.0))
-
-ω1 = ω[i1]
-ue = Xefie[:,i1] / fouriertransform(gaussian)(ω1)
-
-fcr, geo = facecurrents(ue, X)
-Plotly.plot(patch(geo, norm.(fcr)))
-
diff --git a/src/BEAST.jl b/src/BEAST.jl
index 79b1b6ce..d14a5ee2 100644
--- a/src/BEAST.jl
+++ b/src/BEAST.jl
@@ -265,7 +265,6 @@ include("helmholtz3d/timedomain/tdhh3dops.jl")
include("helmholtz3d/timedomain/tdhh3dexc.jl")
include("helmholtz3d/timedomain/tdhh3dpp.jl")
-include("maxwell/timedomain/acustictdops.jl")#support for acusticsinglelayer
include("maxwell/timedomain/mwtdops.jl")
include("maxwell/timedomain/mwtdexc.jl")
include("maxwell/timedomain/tdfarfield.jl")
diff --git a/src/maxwell/timedomain/acustictdops.jl b/src/maxwell/timedomain/acustictdops.jl
deleted file mode 100644
index e25c5078..00000000
--- a/src/maxwell/timedomain/acustictdops.jl
+++ /dev/null
@@ -1,319 +0,0 @@
-mutable struct AcusticSingleLayerTDIO{T} <: RetardedPotential{T}
- "speed of sound in medium"
- speed_of_light::T #for compatibility with the rest of the package I left light instead of sound
- "weight"
- weight::T
- "number of temporal differentiations"
- diffs::Int
-end
-
-function Base.:*(a::Number, op::AcusticSingleLayerTDIO)
- @info "scalar product a * op (acusticSL)"
- AcusticSingleLayerTDIO(
- op.speed_of_light,
- a * op.weight,
- op.diffs)
-end
-
-AcusticSingleLayerTDIO(;speedofsound) = AcusticSingleLayerTDIO(speedofsound, one(speedofsound), 0)
-
-module TDAcustic3D
-import ...BEAST
-
-function acusticsinglelayer(;speedofsound, numdiffs)
- @assert numdiffs >= 0
- #controllare con Kristof che tipo di schema temporale viene utilizzato e come si con integrali rispetto al tempo extra
- #return BEAST.integrate(BEAST.MWSingleLayerTDIO(speedoflight,-1/speedoflight,-speedoflight,2,0))
- return BEAST.AcusticSingleLayerTDIO(speedofsound,one(speedofsound),numdiffs)
-end
-
-end #of the module
-
-export TDAcustic3D
-
-defaultquadstrat(::AcusticSingleLayerTDIO, tfs, bfs) = AllAnalyticalQStrat(1) #nothing
-#nothing goes in hybrid qr, allanalytical goes in zuccottirule
-
-
-function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- testels, trialels, timeels, quadstrat::AllAnalyticalQStrat)
- return nothing
-end
-
-quadrule(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- p, testel, q, trialel, r, timeel, qd, ::AllAnalyticalQStrat) = ZuccottiRule(1.0)
-
-
-
-function quaddata(operator::AcusticSingleLayerTDIO,
- test_local_space, trial_local_space, time_local_space,
- test_element, trial_element, time_element, quadstrat::Nothing)
-
- dmax = numfunctions(time_local_space)-1
- bn = binomial.((0:dmax),(0:dmax)')
-
- V = eltype(test_element[1].vertices)
- ws = WiltonInts84.workspace(V)
- order = 4
- @show order
- quadpoints(test_local_space, test_element, (order,)), bn, ws
-
- end
-
-
- # See: ?BEAST.quadrule for help
- function quadrule(operator::AcusticSingleLayerTDIO,
- test_local_space, trial_local_space, time_local_space,
- p, test_element, q, trial_element, r, time_element,
- quad_data, quadstrat::Nothing)
-
- # WiltonInts84Strat(quad_data[1,p])
- qd = quad_data
- HybridZuccottiWiltonStrat(qd[1][1,p],qd[2],qd[3])
-
- end
-
-
- function innerintegrals!(zlocal, operator::AcusticSingleLayerTDIO,
- test_point,
- test_local_space, trial_local_space, time_local_space,
- test_element, trial_element, time_element,
- quad_rule::HybridZuccottiWiltonStrat, quad_weight)
-
- # error("Here!!!")
-
- dx = quad_weight
- x = cartesian(test_point)
- # n = normal(test_point)
-
- # a = trial_element[1]
- # ξ = x - dot(x -a, n) * n
-
- r = time_element[1]
- R = time_element[2]
- @assert r < R
-
- N = max(degree(time_local_space), 1)
- ∫G, ∫vG, ∫∇G = WiltonInts84.wiltonints(
- trial_element[1],
- trial_element[2],
- trial_element[3],
- x, r, R, Val{2}, quad_rule.workspace)
-
- a = dx / (4*pi)
- #D = operator.num_diffs
- D=0
- @assert D == 0
- @assert numfunctions(test_local_space) == 1
- @assert numfunctions(trial_local_space) == 1
-
- @inline function tmRoR_sl(d, iG)
- sgn = isodd(d) ? -1 : 1
- r = sgn * iG[d+2]
- end
-
- # bns = quad_rule.binomials
-
- @assert D == 0
- for k in 1 : numfunctions(time_local_space)
- d = k - 1
- d < D && continue
- q = reduce(*, d-D+1:d ,init=1)
- zlocal[1,1,k] += a * q * tmRoR_sl(d-D, ∫G)
- end # k
- end
-
-
-
-qpclineline=10^8#quasiparallel case controller
-qpclinetriang=10^10
-qpctriangtriang=10^13
-
-const qpc=[qpclineline, qpclinetriang, qpctriangtriang]
-
-function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::ZuccottiRule) where T
-
- t1=ι[1]
- t2=ι[2]
-
-
- @assert t2 > t1
-
- sos = op.speed_of_light
-
- a1index,a2index,a3index,b1index,b2index,b3index=0,0,0,0,0,0
- Tt = eltype(eltype(τ.vertices))
- hits = 0
- dtol = 1.0e3 * eps(Tt)
- dmin2 = floatmax(Tt)
- for t in 1:3
- for s in 1:3
- d2 = LinearAlgebra.norm_sqr(τ[t]-σ[s])
- dmin2 = min(dmin2, d2)
- if (d2 < dtol)
- hits+=1
- if hits==1
- a1index =t
- b1index = s
- elseif hits==2
- a2index=t
- b2index=s
- end
- end
- end
- end
-
- if hits==3
- #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nt1=",t1,"\nt2=",t2)
-
- z[1,1,1]+=(TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2))/(4*π)
- elseif hits==2
- #pay attention with double layer and index permutation or in whatever case has nx
- if mod1(a1index +1,3)==a2index
- a3index=mod1(a1index-1,3)
- else
- a3index=mod1(a1index+1,3)
- end
-
- if mod1(b1index+1,3)==b2index
- b3index=mod1(b1index-1,3)
- else
- b3index=mod1(b1index+1,3)
- end
- #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
-
-
- z[1,1,1]+=(TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
- elseif hits==1
- a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
- b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
- print("\na1=",τ[a1index],"\na2=",τ[a2index],"\na3=",τ[a3index],"\nb1=",σ[b1index],"\nb2=",σ[b2index],"\nb3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
- z[1,1,1]+=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
-
- else
- #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
- z[1,1,1]+=(inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2,qpc))/(4*π)
- end
- #for the moment sos=1 but I will correct this
-end
-
-
-function momintegrals!(z, op::AcusticSingleLayerTDIO, g::LagrangeRefSpace{T,0,3}, f::LagrangeRefSpace{T,0,3}, t::MonomialBasis{T,0,1}, τ, σ, ι, qr::HybridZuccottiWiltonStrat) where T
-
- t1=ι[1]
- t2=ι[2]
-
-
- @assert t2 > t1
-
- sos = op.speed_of_light
-
- a1index,a2index,a3index,b1index,b2index,b3index=0,0,0,0,0,0
- Tt = eltype(eltype(τ.vertices))
- hits = 0
- dtol = 1.0e3 * eps(Tt)
- dmin2 = floatmax(Tt)
- for t in 1:3
- for s in 1:3
- d2 = LinearAlgebra.norm_sqr(τ[t]-σ[s])
- dmin2 = min(dmin2, d2)
- if (d2 < dtol)
- hits+=1
- if hits==1
- a1index =t
- b1index = s
- elseif hits==2
- a2index=t
- b2index=s
- end
- end
- end
- end
-
- if hits==3
- #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nt1=",t1,"\nt2=",t2)
- entry=(TimeDomainBEMInt.intcoinctriangles(τ[1],τ[2],τ[3],t1,t2))/(4*π)
- if norm(entry)>10.0 || entry<-0.0001
- println("quadrature ",τ[1]," ",τ[2]," ",τ[3]," ", t1," ",t2," ",entry)
- XW = qr.outer_quad_points
- for p in 1 : length(XW)
- x = XW[p].point
- w = XW[p].weight
- innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
- end
- else
- z[1,1,1]+=max(entry,0.0)
- end
- elseif hits==2
- #= XW = qr.outer_quad_points
- for p in 1 : length(XW)
- x = XW[p].point
- w = XW[p].weight
- innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
- end=#
- #pay attention with double layer and index permutation or whatever case has n cross product
- if mod1(a1index +1,3)==a2index
- a3index=mod1(a1index-1,3)
- else
- a3index=mod1(a1index+1,3)
- end
-
- if mod1(b1index+1,3)==b2index
- b3index=mod1(b1index-1,3)
- else
- b3index=mod1(b1index+1,3)
- end
- #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
-
- entry=(TimeDomainBEMInt.inttriangletriangleadjacent(τ[a1index],τ[a2index],τ[a3index],σ[b1index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
- if norm(entry)>10.0 || entry<-0.0001
- println("quadrature ",τ[a1index]," ",τ[a2index]," ",τ[a3index]," ",σ[b1index]," ",σ[b2index]," ",σ[b3index], t1," ",t2," ",entry)
- XW = qr.outer_quad_points
- for p in 1 : length(XW)
- x = XW[p].point
- w = XW[p].weight
- innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
- end
- else
- z[1,1,1]+=max(0.0,entry)
- end
- elseif hits==1
- #=XW = qr.outer_quad_points
- for p in 1 : length(XW)
- x = XW[p].point
- w = XW[p].weight
- innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
- end=#
- a2index,a3index=mod1(a1index+1,3),mod1(a1index+2,3)
- b2index,b3index=mod1(b1index+1,3),mod1(b1index+2,3)
- #print("\np1=",τ[a1index],"\np2=",τ[a2index],"\np3=",τ[a3index],"\nv1=",σ[b1index],"\nv2=",σ[b2index],"\nv3=",σ[b3index],"\nt1=",t1,"\nt2=",t2)
- entry=(TimeDomainBEMInt.intcommonvertex(τ[a1index],τ[a2index],τ[a3index],σ[b2index],σ[b3index],t1,t2,qpc))/(4*π)
- if norm(entry)>10.0 || entry<-0.0001
- println("quadrature ",τ[a1index]," ",τ[a2index]," ",τ[a3index]," ",σ[b2index]," ",σ[b3index], t1," ",t2," ",entry, " ",)
-
- XW = qr.outer_quad_points
- for p in 1 : length(XW)
- x = XW[p].point
- w = XW[p].weight
- innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
- end
- else
- z[1,1,1]+=max(0.0,entry)
- end
- else
- entry=0.1#(inttriangletriangle(τ[1],τ[2],τ[3],σ[1],σ[2],σ[3],t1,t2,qpc))/(4*π)
- #print("\np1=",τ[1],"\np2=",τ[2],"\np3=",τ[3],"\nv1=",σ[1],"\nv2=",σ[2],"\nv3=",σ[3],"\nt1=",t1,"\nt2=",t2)
- if norm(entry)>0.0 || entry<-0.0001
- XW = qr.outer_quad_points
- for p in 1 : length(XW)
- x = XW[p].point
- w = XW[p].weight
- innerintegrals!(z, op, x, g, f, t, τ, σ, ι, qr, w)
- end
-
- else
- z[1,1,1]+=max(0.0,entry)
- end
- end
- #for the moment sos=1 but I will correct this
-end
\ No newline at end of file
diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl
index 632da895..3187a55b 100644
--- a/src/quadrature/quadstrats.jl
+++ b/src/quadrature/quadstrats.jl
@@ -30,9 +30,7 @@ struct OuterNumInnerAnalyticQStrat{R}
outer_rule::R
end
-struct AllAnalyticalQStrat{R}
- rule::R
-end
+
defaultquadstrat(op, tfs, bfs) = defaultquadstrat(op, refspace(tfs), refspace(bfs))
macro defaultquadstrat(dop, body)
diff --git a/src/timedomain/analyticalints.jl b/src/timedomain/analyticalints.jl
deleted file mode 100644
index 97e4359a..00000000
--- a/src/timedomain/analyticalints.jl
+++ /dev/null
@@ -1,232 +0,0 @@
-function minmax1d(vertex,edge)
- T = eltype(τ[1])
- m = norm(vertex-edge[1])
- M = m
- s=edge[1]-edge[2]
- s/=norm(s)
- ev1=edge[1]-vertex
- x0=(edge[1]-dot(ev1,s)*s)
- a=(edge[2]-x0)*s
- b=(edge[1]-x0)*s
- if a<=0 && b>=0
- m=norm(vertex-x0)
- abs(a)<abs(b) && (M=norm(vertex-edge[2]))
-
- else
- for j in 1:length(edge)
- q = edge[j]
- d = norm(vertex-q)
- d < m && (m=d)
- d > M && (M=d)
- end
- end
- return m, M
-end
-
-function rings1d(τ, σ, ΔR)
- m, M = minmaxdist(τ, σ)
- r0 = floor(Int, m/ΔR) + 1
- r1 = ceil(Int, M/ΔR+1)
- r0 : r1
-end
-
-
-
-function quaddata(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- testels::Vector{Simplex{3,0,3,1,T}}, trialels::Vector{Simplex{3,1,1,2,T}}, timeels, quadstrat::AllAnalyticalQStrat) where T
-
- dimU=dimension(testels)
- dimV=dimension(trialels)
- #rigerenerare delta R
-#quaddata 1D
- #if dimU+dimV==1
- #testelsboundary=skeleton(testels,dimU-1)
- # trialelsboundary=skeleton(trialels,dimV-1)
-
- numnodes=length(testels)
- numedges=length(trialels)
-
- datavertexedge=Array{TimeDomainBEMInt.edgevertexgeo{T,P}, 2}(undef, numnodes, numedges)
- rings=Array{UnitRange{Int},2}(undef, numnodes, numedges)
- datarings=Array{Vector{Tuple{Int,Vector}},2}(undef,numnodes,numedges)#il type va bene
-
- #fill datarings with zeross !
-
- for p in 1:numnodes
- τ = chart(testels,p)#testels[p]
- for q in 1:numedges
- σ = chart(trialels,q)
- edgevertgeo=TimeDomainBEMInt.edgevertexinteraction(τ,σ[1],σ[2])
- datavertexedge[p,q]=edgevertgeo
- a,b=edgevertegeo.extint0[1],edgevertegeo.extint0[2]
- rngs=rings1d(τ,σ,ΔR)
- rings[p,q]=rngs
- datarings[p,q]=[0,[0.0,0.0]]
- for r in rngs
- r > numfunctions(timebasisfunction) && continue #serve?
- ι = ring(r,ΔR)
-
- # compute interactions between reference shape functions
- #fill!(z, 0)
- rp=τ #se e un simplex ok se no va messo chart(τ,1).vertices credo
- t2=ι[2]#needs a check
- extint=TimeDomainBEMInt.edgevertexinteraction(t2,edgevertgeo)
- push!(datarings[p,q],extint)
-
- # qr = quadrule(op, U, V, W, p, τ, q, σ, r, ι, qd, quadstrat)
- #momintegrals!(z, op, U, V, W, τ, σ, ι, qr)
- end
- end
- end
-
- return datavertexedge
- #else
- # return "devo ancora scrivere"
- #end
-end
-
-function quaddata2Dee(op::AcusticSingleLayerTDIO, testrefs, trialrefs, timerefs,
- testels::Vector{Simplex{3,1,1,2,T}}, trialels::Vector{Simplex{3,1,1,2,T}}, timeels, quadstrat::AllAnalyticalQStrat)
- nnodes=length(nodes)
-
-
- totrings=Array{UnitRange{Int},2}(undef, numedges, numedges)
- datavalues=Array{Vector{Tuple{Int,Vector}},2}(undef,numedges,numedges)
- cnnct=connectivity(edges,nodes)
-
- for p in 1:numdges
- for q in 1:numedges
- edge1=chart(edges,p)
- edge2=chart(edges,q)
-
- vertind1=cnnct[1:nnodes,p].nzind
- vertsgn1=cnnct[1:nnodes,p].nzval
- vertind2=cnnct[1:nnodes,q].nzind
- vertsgn2=cnnct[1:nnodes,q].nzval
-
- if vertsgn1[1]==1
- a1,a2=edge1[1],edge1[2]
- else
- a2,a1=edge1[1],edge1[2]
- end
-
- if vertsgn2[1]==1
- b1,b2=edge2[1],edge2[2]
- else
- b2,b1=edge2[1],edge2[2]
- end
-
- geo1,rings1,datarings1=edgevertexgeo[vertind1[1],q],rings[vertind1[1],q],datarings[vertind1[1],q]
- geo2,rings2,datarings2=edgevertexgeo[vertind1[2],q],rings[vertind1[2],q],datarings[vertind1[2],q]
- geo3,rings3,datarings3=edgevertexgeo[vertind2[1],p],rings[vertind2[1],p],datarings[vertind2[1],p]
- geo4,rings4,datarings4=edgevertexgeo[vertind2[2],p],rings[vertind2[2],p],datarings[vertind2[2],p]
-
- geo=[geo1,geo2,geo3,geo4]
- rings=[rings1,rings2,rings3,rings4]
- datarings=[datarings1,datarings2,datarings3,datarings4]
-
- totrings[p,q],datavalues[p,q]=intlinelineglobal(a1,a2,b1,b2,geo,rings,datarings,[10^6,10^6,10^6],Val{0})
- end
- end
- return totrings,datavalues
-end
-
-
-function intlinelineglobal(a1,a2,b1,b2,geo,rings,datarings,parcontrol,UB::Type{Val{N}}) where N
-
- #nedges=length(edges)
-
-
- #vertices=[a1,a2,a1′,a2′]
- vertices=[a1,a2,b1,b2]
- l12=norm(vertices[1]-vertices[2])
- l12′=norm(vertices[3]-vertices[4])
-
-
- #geo1=edgevertexinteraction(a1,a1′,a2′)
- #geo2=edgevertexinteraction(a2,a1′,a2′)
- #geo3=edgevertexinteraction(a1′,a1,a2)
- #geo4=edgevertexinteraction(a2′,a1,a2)
-
- #datatime=Array{Tuple}(undef,2,4)?
- I = maketuple(eltype(a1), UB)
- K = maketuple(typeof(a1), UB)
-
- x=geo[3].tangent
-
- #z=cross(a12′,x)
- #J=norm(z)
- #if J blablabla
- #z /= J
- #h=dot(r22,z)
- hdir=cross(geo[1].tangent,x)
- n=hdir/norm(hdir)
- sgnn=[+1,-1,-1,+1]
- h=dot(a2-b2,n)
- sgnh=[+1,-1,+1,-1]
- angletot=0.0
- dminv=Vector{eltype(edge1[1])}(undef, 4)
- ξ=Vector{typeof(edge1[1])}(undef, 4)
- for j in 1:4
- dminv[j]=geo[j].dmin
- dmaxv[j]=geo[j].dmax
- v=vertices[j]
- ξ[j]=v-n*h*sgnh[j]*sgnn[j]
- angletot+=anglecontribution(ξ[j],sgnn[j]*n,geo[j])
- end
- if abs(angletot-2π)<100*eps(eltype(edge1[1]))
- dmin=abs(h)
- else
- dmin=min(dminv[1],dminv[2],dminv[3],dminv[4])
- end
-
- dmax=max(dmaxv[1],dmaxv[2],dmaxv[3],dmaxv[4])
-
- r0 = floor(Int, dmin/ΔR) + 1
- r1 = ceil(Int, dmax/ΔR+1) #recuperare deltaR
- ringtot = r0 : r1
-
- allint=Vector{typeof((I,K))}(undef,r1-r0+2)
- fill!(allint,(I,K))
- if norm(hdir) < (parcontrol[1])*eps(typeof(temp1))
- I=intparallelsegment(a1,a2,b1,b2,temp1,temp2)[1] #attenzione qui non compatibile con quello che stiamo scrivendo
- else
- n=hdir/norm(hdir)
- sgnn=[+1,-1,-1,+1]
- h=dot(a2-a2′,n)
- sgnh=[+1,-1,+1,-1]
- for j in 1:4
- for i in ringtot[1]:(rings[j][1]-1)
-
- P,Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],0,[0,0],i*ΔR,UB)
- allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
- allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
- P,Q = arcsegcontribution(v,ξ[j],-sgnn[j]*n,sgnh[j]*h,geo[j],0,[0,0],(i-1)*ΔR,UB)
- allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
- allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
- end
- for i in rings[j]
-
- #shall I put some check like i*deltaR > h
- P, Q = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-rings[j][1]+2][1],datarings[j][i-rings[j][1]+2][2],i*ΔR,UB) #salviamo
- allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
- allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
- P, Q = arcsegcontribution(v,ξ[j],-sgnn[j]*n,sgnh[j]*h,geo[j],datarings[j][i-rings[j][1]+1][1],datarings[j][i-rings[j][1]+1][2],(i-1)*ΔR,UB)
- allint[i-ringtot[1]+2][1] = add(allint[i-ringtot[1]+2][1],P)
- allint[i-ringtot[1]+2][2] = add(allint[i-ringtot[1]+2][2],Q)
- end #probabilmente va bene cosi anche con ceil(int,frac) invece di ceil(int,frac+1)
- allint[i-ringtot[1]+2][1]=multiply(allint[i-ringtot[1]+2][1],1/(l12′*l12*norm(hdir)))
- allint[i-ringtot[1]+2][2]=multiply(allint[i-ringtot[1]+2][2],1/(l12′*l12*norm(hdir)))
- #= for i in (relrings[j][2]+1):ringtot[j][2]
-
- #shall I put some check like i*deltaR > h
- saveP[i],saveQ[i] = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],1,[a,b],i*ΔR,UB) #dadefinire a e b
- save,and,subtract = arcsegcontribution(v,ξ[j],sgnn[j]*n,sgnh[j]*h,geo[j],1,[a,b],(i-1)*ΔR,UB)
- end =# #sembra che non serva a causa di ceil(int,frac+1)!
- end
- #I+=(1/abs(r12′[2]*r12[1]))*(3*(temp2^2-temp1^2)*d[1]-2*(temp2^3-temp1^3)*d[2])
- end
-
- return ringtot,allint #missing buidgrad since it is not yet adapted for int line line
-end
-
diff --git a/src/timedomain/tdintegralop.jl b/src/timedomain/tdintegralop.jl
index 9834e5bc..8ec3cd44 100644
--- a/src/timedomain/tdintegralop.jl
+++ b/src/timedomain/tdintegralop.jl
@@ -1,5 +1,4 @@
using WiltonInts84
-using TimeDomainBEMInt
abstract type AbstractSpaceTimeOperator end
abstract type SpaceTimeOperator <: AbstractSpaceTimeOperator end # atomic operator
@@ -315,15 +314,6 @@ struct WiltonInts84Strat{T,V,W}
workspace::W
end
-struct HybridZuccottiWiltonStrat{T,V,W}
- outer_quad_points::T
- binomials::V
- workspace::W
-end
-
-struct ZuccottiRule{T}
- weight::T #Allanalyticalformula
-end
function momintegrals!(z, op, g, f, T, τ, σ, ι, qr::WiltonInts84Strat)
From 7537d702461c9c96c95884e846dcc60a6cd66345 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 9 Oct 2024 16:38:45 +0200
Subject: [PATCH 15/24] Update quadstrats.jl
---
src/quadrature/quadstrats.jl | 1 -
1 file changed, 1 deletion(-)
diff --git a/src/quadrature/quadstrats.jl b/src/quadrature/quadstrats.jl
index 3187a55b..534f9a46 100644
--- a/src/quadrature/quadstrats.jl
+++ b/src/quadrature/quadstrats.jl
@@ -31,7 +31,6 @@ struct OuterNumInnerAnalyticQStrat{R}
end
-
defaultquadstrat(op, tfs, bfs) = defaultquadstrat(op, refspace(tfs), refspace(bfs))
macro defaultquadstrat(dop, body)
@assert dop.head == :tuple
From 95789fa4a7bdd232d7a50f4ecb53ffb5148a6cb9 Mon Sep 17 00:00:00 2001
From: bmergl <bmergl@outlook.de>
Date: Thu, 24 Oct 2024 18:39:40 +0200
Subject: [PATCH 16/24] VIE bugfix
Fix for VIE momintegrals! and grideval.
(This was needed due to changes in momintegrals! and numfunctions.)
Activates the VIE unit test.
---
src/bases/local/laglocal.jl | 5 +++--
src/postproc/segcurrents.jl | 12 ++++++++----
src/volumeintegral/sauterschwab_ints.jl | 19 +++++++++++--------
test/runtests.jl | 2 ++
test/test_hh_lsvie.jl | 18 +++++++++---------
5 files changed, 33 insertions(+), 23 deletions(-)
diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl
index 21355e0f..0adea811 100644
--- a/src/bases/local/laglocal.jl
+++ b/src/bases/local/laglocal.jl
@@ -9,8 +9,9 @@ numfunctions(s::LagrangeRefSpace{T,1,3}, ch::CompScienceMeshes.ReferenceSimplex{
numfunctions(s::LagrangeRefSpace{T,2,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 6
numfunctions(s::LagrangeRefSpace{T,Dg}, ch::CompScienceMeshes.ReferenceSimplex{D}) where {T,Dg,D} = binomial(D+Dg,Dg)
-valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
- SVector{numfunctions(ref), Tuple{T,T}}
+# valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} =
+# SVector{numfunctions(ref), Tuple{T,T}}
+valuetype(ref::LagrangeRefSpace{T}, charttype) where {T} = T
# Evaluate constant lagrange elements on anything
(ϕ::LagrangeRefSpace{T,0})(tp) where {T} = SVector(((value=one(T), derivative=zero(T)),))
diff --git a/src/postproc/segcurrents.jl b/src/postproc/segcurrents.jl
index a02e37c9..844f2eec 100644
--- a/src/postproc/segcurrents.jl
+++ b/src/postproc/segcurrents.jl
@@ -29,21 +29,25 @@ function grideval(points, coeffs, basis; type=nothing)
V = valuetype(refs, eltype(charts))
T = promote_type(eltype(coeffs), eltype(V))
- P = similar_type(V, T)
+ if !(V <: SVector)
+ P = T
+ else
+ P = similar_type(V, T)
+ end
- type != nothing && (P = type)
+ type !== nothing && (P = type)
values = zeros(P, size(points))
chart_tree = BEAST.octree(charts)
for (j,point) in enumerate(points)
i = CompScienceMeshes.findchart(charts, chart_tree, point)
- if i != nothing
+ if i !== nothing
# @show i
chart = charts[i]
u = carttobary(chart, point)
vals = refs(neighborhood(chart,u))
- for r in 1 : numfunctions(refs)
+ for r in 1 : numfunctions(refs, domain(chart))
for (m,w) in ad[i, r]
values[j] += w * coeffs[m] * vals[r][1]
end
diff --git a/src/volumeintegral/sauterschwab_ints.jl b/src/volumeintegral/sauterschwab_ints.jl
index 68385e7a..6b8fc46d 100644
--- a/src/volumeintegral/sauterschwab_ints.jl
+++ b/src/volumeintegral/sauterschwab_ints.jl
@@ -66,7 +66,7 @@ function reorder_dof(space::RTRefSpace,I)
return SVector(K),SVector{3,Int64}(1,1,1)
end
-function reorder_dof(space::LagrangeRefSpace{Float64,0,3,1},I)
+function reorder_dof(space::LagrangeRefSpace{T,0,3,1},I) where T
return SVector{1,Int64}(1),SVector{1,Int64}(1)
end
@@ -102,16 +102,19 @@ function reorder_dof(space::LagrangeRefSpace{T,1,4,4},I) where T
end
function momintegrals!(out, op::VIEOperator,
- test_local_space::RefSpace, test_ptr, test_tetrahedron_element,
- trial_local_space::RefSpace, trial_ptr, trial_tetrahedron_element,
+ test_functions::Space, test_ptr, test_tetrahedron_element,
+ trial_functions::Space, trial_ptr, trial_tetrahedron_element,
strat::SauterSchwab3DStrategy)
+ local_test_space = refspace(test_functions)
+ local_trial_space = refspace(trial_functions)
+
#Find permutation of vertices to match location of singularity to SauterSchwab
J, I= SauterSchwab3D.reorder(strat.sing)
#Get permutation and rel. orientatio of DoFs
- K,O1 = reorder_dof(test_local_space, I)
- L,O2 = reorder_dof(trial_local_space, J)
+ K,O1 = reorder_dof(local_test_space, I)
+ L,O2 = reorder_dof(local_trial_space, J)
#Apply permuation to elements
if length(I) == 4
@@ -146,7 +149,7 @@ function momintegrals!(out, op::VIEOperator,
#Define integral (returns a function that only needs barycentric coordinates)
igd = VIEIntegrand(test_tetrahedron_element, trial_tetrahedron_element,
- op, test_local_space, trial_local_space)
+ op, local_test_space, local_trial_space)
#Evaluate integral
Q = SauterSchwab3D.sauterschwab_parameterized(igd, strat)
@@ -161,8 +164,8 @@ function momintegrals!(out, op::VIEOperator,
end
function momintegrals!(z, biop::VIEOperator,
- tshs, tptr, tcell,
- bshs, bptr, bcell,
+ test_functions::Space, tptr, tcell,
+ trial_functions::Space, bptr, bcell,
strat::DoubleQuadRule)
# memory allocation here is a result from the type instability on strat
diff --git a/test/runtests.jl b/test/runtests.jl
index bbcee67f..f99b04bf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -79,6 +79,8 @@ include("test_variational.jl")
include("test_handlers.jl")
include("test_gridfunction.jl")
+include("test_hh_lsvie.jl")
+
using TestItemRunner
@run_package_tests
diff --git a/test/test_hh_lsvie.jl b/test/test_hh_lsvie.jl
index 899c55b1..cc16cb35 100644
--- a/test/test_hh_lsvie.jl
+++ b/test/test_hh_lsvie.jl
@@ -11,12 +11,12 @@ using Test
@testset "Lippmann Schwinger Volume Integral Equation" begin
# Environment
- ε1 = 1.0*ε0
- μ1 = μ0
+ ε1 = 1.0*SphericalScattering.ε0
+ μ1 = SphericalScattering.μ0
# Dielectic Sphere
- ε2 = 5.0*ε0
- μ2 = μ0
+ ε2 = 5.0*SphericalScattering.ε0
+ μ2 = SphericalScattering.μ0
r = 1.0
sp = DielectricSphere(; radius = r, filling = Medium(ε2, μ2))
@@ -56,7 +56,7 @@ using Test
# Assembly
- b = assemble(Φ_inc, X)
+ b = real.(assemble(Φ_inc, X))
Z_I = assemble(I, X, X)
@@ -69,8 +69,8 @@ using Test
Z_version2 = Z_I + Z_Y
# MoM solution
- u_version1 = Z_version1 \ Vector(b)
- u_version2 = Z_version2 \ Vector(b)
+ u_version1 = Z_version1 \ b
+ u_version2 = Z_version2 \ b
# Observation points
range_ = range(-1.0*r,stop=1.0*r,length=14)
@@ -84,8 +84,8 @@ using Test
Φ = field(sp, ex, ScalarPotential(points_sp))
# MoM solution inside the dielectric sphere
- Φ_MoM_version1 = BEAST.grideval(points_sp, u_version1, X, type=Float64)
- Φ_MoM_version2 = BEAST.grideval(points_sp, u_version2, X, type=Float64)
+ Φ_MoM_version1 = BEAST.grideval(points_sp, u_version1, X)
+ Φ_MoM_version2 = BEAST.grideval(points_sp, u_version2, X)
From 4d485d145644e8de9acfa1ca244ad00294d510b1 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Tue, 12 Nov 2024 15:04:36 +0100
Subject: [PATCH 17/24] Update builddual.jl
---
examples/builddual.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/examples/builddual.jl b/examples/builddual.jl
index 84308adc..cba47e7c 100644
--- a/examples/builddual.jl
+++ b/examples/builddual.jl
@@ -3,7 +3,7 @@ using BEAST
Faces = meshsphere(1.0, 0.35)
Edges = skeleton(Faces,1)
-
+#ciao
faces = barycentric_refinement(Faces)
edges = skeleton(faces,1)
# verts = skeleton(faces,0)
From d234b221cf423a3edab29effd692615bac8769ad Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Tue, 12 Nov 2024 15:06:23 +0100
Subject: [PATCH 18/24] Update builddual.jl
---
examples/builddual.jl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/examples/builddual.jl b/examples/builddual.jl
index cba47e7c..84308adc 100644
--- a/examples/builddual.jl
+++ b/examples/builddual.jl
@@ -3,7 +3,7 @@ using BEAST
Faces = meshsphere(1.0, 0.35)
Edges = skeleton(Faces,1)
-#ciao
+
faces = barycentric_refinement(Faces)
edges = skeleton(faces,1)
# verts = skeleton(faces,0)
From f387b8a7f73f7819e905e76aca2fbae0268d2112 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Tue, 12 Nov 2024 15:25:49 +0100
Subject: [PATCH 19/24] Update Project.toml
---
Project.toml | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/Project.toml b/Project.toml
index a898b184..1b31d623 100644
--- a/Project.toml
+++ b/Project.toml
@@ -47,11 +47,11 @@ FillArrays = "0.11, 0.12, 0.13, 1"
Infiltrator = "1.8.2"
IterativeSolvers = "0.9"
LiftedMaps = "0.5.1"
-LinearMaps = "3.7 - 3.9, 3"
+LinearMaps = "3.7 - 3.9, 3.11.2"
NestedUnitRanges = "0.2"
Requires = "1"
SauterSchwab3D = "0.1"
-SauterSchwabQuadrature = "2.3.0"
+SauterSchwabQuadrature = "2.4.0"
SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1, 2"
StaticArrays = "0.8.3, 0.9, 0.10, 0.11, 0.12, 1"
TestItems = "0.1.1, 1"
From 4fa214fea3a5d521ed5d92a294026dd5b56e0d73 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 13 Nov 2024 16:55:58 +0100
Subject: [PATCH 20/24] Update ncrossbdmlocal.jl
---
src/bases/local/ncrossbdmlocal.jl | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/src/bases/local/ncrossbdmlocal.jl b/src/bases/local/ncrossbdmlocal.jl
index 72bd24e5..8ec3dced 100644
--- a/src/bases/local/ncrossbdmlocal.jl
+++ b/src/bases/local/ncrossbdmlocal.jl
@@ -19,6 +19,8 @@ function (f::NCrossBDMRefSpace{T})(p) where T
(value= n × (v*tv)/j, curl=d),]
end
+numfunctions(x::NCrossBDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
+#=
const _vert_perms_ncrossbdm = [
(1,2,3),
(2,3,1),
@@ -35,7 +37,7 @@ const _dof_perms_ncrossbdm = [
(2,1,6,5,4,3),
(6,5,4,3,2,1),
]
-
+=#
function dof_permutation(::NCrossBDMRefSpace, vert_permutation)
i = findfirst(==(tuple(vert_permutation...)), _vert_perms_ncrossbdm)
return _dof_perms_ncrossbdm[i]
From 809e480beecc158035d58a8d5fc8ada5444a6264 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 13 Nov 2024 17:36:20 +0100
Subject: [PATCH 21/24] perms of curl conf spaces not needed
---
src/bases/local/bdmlocal.jl | 2 ++
src/bases/local/ncrossbdmlocal.jl | 22 ----------------------
src/bases/local/ndlocal.jl | 23 -----------------------
3 files changed, 2 insertions(+), 45 deletions(-)
diff --git a/src/bases/local/bdmlocal.jl b/src/bases/local/bdmlocal.jl
index 1c83f5f7..605963d9 100644
--- a/src/bases/local/bdmlocal.jl
+++ b/src/bases/local/bdmlocal.jl
@@ -21,6 +21,8 @@ end
divergence(ref::BDMRefSpace, sh, el) = [Shape(sh.cellid, 1, sh.coeff/(2*volume(el)))]
+numfunctions(x::BDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
+
const _vert_perms_bdm = [
(1,2,3),
(2,3,1),
diff --git a/src/bases/local/ncrossbdmlocal.jl b/src/bases/local/ncrossbdmlocal.jl
index 8ec3dced..cb332499 100644
--- a/src/bases/local/ncrossbdmlocal.jl
+++ b/src/bases/local/ncrossbdmlocal.jl
@@ -20,25 +20,3 @@ function (f::NCrossBDMRefSpace{T})(p) where T
end
numfunctions(x::NCrossBDMRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 6
-#=
-const _vert_perms_ncrossbdm = [
- (1,2,3),
- (2,3,1),
- (3,1,2),
- (2,1,3),
- (1,3,2),
- (3,2,1),
-]
-const _dof_perms_ncrossbdm = [
- (1,2,3,4,5,6),
- (5,6,1,2,3,4),
- (3,4,5,6,1,2),
- (4,3,2,1,6,5),
- (2,1,6,5,4,3),
- (6,5,4,3,2,1),
-]
-=#
-function dof_permutation(::NCrossBDMRefSpace, vert_permutation)
- i = findfirst(==(tuple(vert_permutation...)), _vert_perms_ncrossbdm)
- return _dof_perms_ncrossbdm[i]
-end
\ No newline at end of file
diff --git a/src/bases/local/ndlocal.jl b/src/bases/local/ndlocal.jl
index 93a31471..a4b3bd17 100644
--- a/src/bases/local/ndlocal.jl
+++ b/src/bases/local/ndlocal.jl
@@ -68,26 +68,3 @@ function restrict(ϕ::NDRefSpace{T}, dom1, dom2) where T
return Q
end
-
-const _vert_perms_nd = [
- (1,2,3),
- (2,3,1),
- (3,1,2),
- (2,1,3),
- (1,3,2),
- (3,2,1),
-]
-const _dof_perms_nd = [
- (1,2,3),
- (3,1,2),
- (2,3,1),
- (2,1,3),
- (1,3,2),
- (3,2,1),
-]
-
-function dof_permutation(::NDRefSpace, vert_permutation)
- i = findfirst(==(tuple(vert_permutation...)), _vert_perms_nd)
- @assert i != nothing
- return _dof_perms_nd[i]
-end
\ No newline at end of file
From 105923207aab9140260b9a6188d50e40356932b0 Mon Sep 17 00:00:00 2001
From: azuccott <alessandro.zuccotti@ugent.be>
Date: Wed, 13 Nov 2024 17:51:30 +0100
Subject: [PATCH 22/24] space correction
---
src/maxwell/mwops.jl | 2 --
src/timedomain/tdintegralop.jl | 2 --
2 files changed, 4 deletions(-)
diff --git a/src/maxwell/mwops.jl b/src/maxwell/mwops.jl
index 3e8429b4..dbac76e5 100644
--- a/src/maxwell/mwops.jl
+++ b/src/maxwell/mwops.jl
@@ -82,8 +82,6 @@ singularpart(op::MWDoubleLayer3D) = MWDoubleLayer3DSng(op.alpha, op.gamma)
-
-
# function quadrule(op::MaxwellOperator3D, g::BDMRefSpace, f::BDMRefSpace, i, τ, j, σ, qd,
# qs::DoubleNumWiltonSauterQStrat)
diff --git a/src/timedomain/tdintegralop.jl b/src/timedomain/tdintegralop.jl
index a0c1bee9..73bcb201 100644
--- a/src/timedomain/tdintegralop.jl
+++ b/src/timedomain/tdintegralop.jl
@@ -328,5 +328,3 @@ function momintegrals!(z, op, g, f, T, τ, σ, ι, qr::WiltonInts84Strat)
end
end
-
-
From 263926d826fdff3475248042411ba791907a611e Mon Sep 17 00:00:00 2001
From: Kristof Cools <krcools@gmail.com>
Date: Fri, 20 Dec 2024 11:07:43 +0100
Subject: [PATCH 23/24] some tests got unborken by update in WiltonInts
---
test/test_tdmwdbl.jl | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/test/test_tdmwdbl.jl b/test/test_tdmwdbl.jl
index 81fb700f..60d6c453 100644
--- a/test/test_tdmwdbl.jl
+++ b/test/test_tdmwdbl.jl
@@ -66,7 +66,7 @@ qrl = BEAST.quadrule(K, refspace(X), refspace(Y), refspace(T2),
1, test_els[1], 1, trial_els[1], 1, time_els[1], qdt, qs)
z = zeros(Float64, 3, 3, 10)
BEAST.innerintegrals!(z, K, x, refspace(X), refspace(Y), refspace(T2), test_els[1], trial_els[1], (0.0, 20.0), qrl, w)
-@test_broken !any(isnan.(z))
+@test !any(isnan.(z))
verts = trial_els[1].vertices
import WiltonInts84
@@ -85,4 +85,4 @@ h = 0.0
using StaticArrays
m = SVector(-0.7071067811865474, 0.7071067811865478, -0.0)
iG, vG = WiltonInts84.segintsg(a, b, p, h, m, Val{2})
-@test_broken !any(isnan.(iG))
+@test !any(isnan.(iG))
From 0f4fa419ac787ba3e1bd03a98c7bdf60d8fd8ea0 Mon Sep 17 00:00:00 2001
From: Kristof Cools <krcools@gmail.com>
Date: Fri, 20 Dec 2024 11:40:23 +0100
Subject: [PATCH 24/24] unbreak tests in tdmwdbl
---
test/test_tdmwdbl.jl | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/test/test_tdmwdbl.jl b/test/test_tdmwdbl.jl
index 81fb700f..60d6c453 100644
--- a/test/test_tdmwdbl.jl
+++ b/test/test_tdmwdbl.jl
@@ -66,7 +66,7 @@ qrl = BEAST.quadrule(K, refspace(X), refspace(Y), refspace(T2),
1, test_els[1], 1, trial_els[1], 1, time_els[1], qdt, qs)
z = zeros(Float64, 3, 3, 10)
BEAST.innerintegrals!(z, K, x, refspace(X), refspace(Y), refspace(T2), test_els[1], trial_els[1], (0.0, 20.0), qrl, w)
-@test_broken !any(isnan.(z))
+@test !any(isnan.(z))
verts = trial_els[1].vertices
import WiltonInts84
@@ -85,4 +85,4 @@ h = 0.0
using StaticArrays
m = SVector(-0.7071067811865474, 0.7071067811865478, -0.0)
iG, vG = WiltonInts84.segintsg(a, b, p, h, m, Val{2})
-@test_broken !any(isnan.(iG))
+@test !any(isnan.(iG))