From b12990ba82b155e0e44a1291ee58b33e4e49e634 Mon Sep 17 00:00:00 2001 From: Kristof Cools Date: Tue, 1 Oct 2024 14:44:10 +0200 Subject: [PATCH] numfunctions now requires both the refspace and refdomain --- src/BEAST.jl | 17 ++-- src/bases/basis.jl | 18 ++-- src/bases/global/globaldofs.jl | 19 +--- src/bases/local/bdm3dlocal.jl | 2 +- src/bases/local/bdmlocal.jl | 2 +- src/bases/local/gwplocal.jl | 76 +++++++++++++++ src/bases/local/laglocal.jl | 20 ++-- src/bases/local/ncrossbdmlocal.jl | 2 +- src/bases/local/nd2local.jl | 4 +- src/bases/local/ndlcclocal.jl | 6 +- src/bases/local/ndlcdlocal.jl | 6 +- src/bases/local/ndlocal.jl | 5 +- src/bases/local/rt2local.jl | 3 +- src/bases/local/rtlocal.jl | 4 +- src/bases/local/rtqlocal.jl | 2 +- src/bases/localbasis.jl | 29 ++++++ src/bases/subdbasis.jl | 2 +- src/bases/timebasis.jl | 4 +- src/excitation.jl | 10 +- src/helmholtz3d/timedomain/tdhh3dops.jl | 24 ++--- src/helmholtz3d/wiltonints.jl | 92 +++++++++++++------ src/integralop.jl | 31 +++++-- src/localop.jl | 37 +++++++- src/maxwell/timedomain/mwtdops.jl | 14 ++- src/maxwell/wiltonints.jl | 14 ++- src/postproc.jl | 4 +- src/quadrature/rules/testrefinestrialqrule.jl | 10 +- src/quadrature/rules/trialrefinestestqrule.jl | 10 +- src/quadrature/sauterschwabints.jl | 5 +- src/timedomain/tdexcitation.jl | 16 +++- src/timedomain/tdintegralop.jl | 7 +- src/utils/lagpolys.jl | 45 +++++++++ test/test_basis.jl | 6 +- test/test_gwp.jl | 12 +++ test/test_rt.jl | 3 +- test/test_ss_nested_meshes.jl | 28 +++--- test/test_tdassembly.jl | 8 +- test/test_tdhhdbl.jl | 6 +- test/test_wiltonints.jl | 2 +- 39 files changed, 449 insertions(+), 156 deletions(-) create mode 100644 src/bases/local/gwplocal.jl create mode 100644 src/bases/localbasis.jl create mode 100644 src/utils/lagpolys.jl create mode 100644 test/test_gwp.jl diff --git a/src/BEAST.jl b/src/BEAST.jl index 10568c72..96d4fc5d 100644 --- a/src/BEAST.jl +++ b/src/BEAST.jl @@ -141,14 +141,9 @@ include("utils/combinatorics.jl") include("utils/linearspace.jl") include("utils/zeromap.jl") include("utils/rank1map.jl") +include("utils/lagpolys.jl") -include("bases/basis.jl") -include("bases/lincomb.jl") -include("bases/trace.jl") -include("bases/restrict.jl") -include("bases/divergence.jl") -include("bases/global/globaldofs.jl") - +include("bases/localbasis.jl") include("bases/local/laglocal.jl") include("bases/local/rtlocal.jl") include("bases/local/rt2local.jl") @@ -160,6 +155,14 @@ include("bases/local/ndlcclocal.jl") include("bases/local/ndlcdlocal.jl") include("bases/local/bdm3dlocal.jl") include("bases/local/rtqlocal.jl") +include("bases/local/gwplocal.jl") + +include("bases/basis.jl") +include("bases/lincomb.jl") +include("bases/trace.jl") +include("bases/restrict.jl") +include("bases/divergence.jl") +include("bases/global/globaldofs.jl") include("bases/lagrange.jl") include("bases/rtspace.jl") diff --git a/src/bases/basis.jl b/src/bases/basis.jl index e2bee21d..e75a0733 100644 --- a/src/bases/basis.jl +++ b/src/bases/basis.jl @@ -1,5 +1,4 @@ -abstract type RefSpace{T,D} end -abstract type DivRefSpace{T,D} <: RefSpace{T,D} end + abstract type AbstractSpace end abstract type Space{T} <: AbstractSpace end @@ -25,12 +24,12 @@ Returns the ReferenceSpace of local shape functions on which the basis is built. function refspace end -""" - numfunctions(r) +# """ +# numfunctions(r) -Return the number of functions in a `Space` or `RefSpace`. -""" -numfunctions(rs::RefSpace{T,D}) where {T,D} = D +# Return the number of functions in a `Space` or `RefSpace`. +# """ +# numfunctions(rs::RefSpace{T,D}) where {T,D} = D """ @@ -208,7 +207,10 @@ function assemblydata(basis::Space; onlyactives=true) num_cells = numcells(geo) num_bfs = numfunctions(basis) - num_refs = numfunctions(refspace(basis)) + + ch = chart(geo, first(geo)) + dom = domain(ch) + num_refs = numfunctions(refspace(basis), dom) # Determine the number of functions defined over a given cell # and per local shape function. diff --git a/src/bases/global/globaldofs.jl b/src/bases/global/globaldofs.jl index be974975..b736bb1c 100644 --- a/src/bases/global/globaldofs.jl +++ b/src/bases/global/globaldofs.jl @@ -20,16 +20,7 @@ function trace(edge_ch, face_ch, localspace::RefSpace) return f end -function _lagrangepolynomial(nodes, i, s, i1=length(nodes)) - r = one(T) - si = nodes[i] - for j in 1:i1 - j == i && continue - sj = nodes[j] - r *= (s - sj) / (si - sj) - end - return r -end + struct _LagrangeGlobalEdgeDoFs order::Int @@ -38,7 +29,7 @@ function numfunctions(dof::_LagrangeGlobalEdgeDoFs) dof.order-1 end function (dof::_LagrangeGlobalEdgeDoFs)(s) T = typeof(s) nodes = range(zero(T), one(T), length=order+1) - [_lagrangepolynomial(nodes, i, s) for i in 2:order] + [_lagpoly(nodes, i, s) for i in 2:order] end struct _LagrangeGlobalFaceDoFs @@ -54,12 +45,12 @@ end # idx = 1 # degree = dof.degree # for i in 0:degree -# prodi = _lagrangepolynomial(nodes, i+1, s1, i) +# prodi = _lagpoly(nodes, i+1, s1, i) # for j in 0:degree # k = degree - i - j # k < 0 && continue -# prodj = _lagrangepolynomial(nodes, j+1, s2, j) -# prodk = _lagrangepolynomial(nodes, k+1, s3, k) +# prodj = _lagpoly(nodes, j+1, s2, j) +# prodk = _lagpoly(nodes, k+1, s3, k) # r[idx] = prodi * prodj * prodk # end end # return r diff --git a/src/bases/local/bdm3dlocal.jl b/src/bases/local/bdm3dlocal.jl index 9b1694ea..02b3b02a 100644 --- a/src/bases/local/bdm3dlocal.jl +++ b/src/bases/local/bdm3dlocal.jl @@ -1,4 +1,4 @@ -struct BDM3DRefSpace{T} <: RefSpace{T,12} end +struct BDM3DRefSpace{T} <: RefSpace{T} end function (f::BDM3DRefSpace)(p) diff --git a/src/bases/local/bdmlocal.jl b/src/bases/local/bdmlocal.jl index a07e42b6..1eb173a4 100644 --- a/src/bases/local/bdmlocal.jl +++ b/src/bases/local/bdmlocal.jl @@ -1,4 +1,4 @@ -struct BDMRefSpace{T} <: RefSpace{T,6} end +struct BDMRefSpace{T} <: RefSpace{T} end function (f::BDMRefSpace)(p) diff --git a/src/bases/local/gwplocal.jl b/src/bases/local/gwplocal.jl new file mode 100644 index 00000000..a0788f4a --- /dev/null +++ b/src/bases/local/gwplocal.jl @@ -0,0 +1,76 @@ +struct GWPCurlRefSpace{T,Degree} <: RefSpace{T} end + +function (ϕ::GWPCurlRefSpace{T,Degree})(p) where {T,Degree} + dom = domain(chart(p)) + u = parametric(p) + vals = ϕ(dom, u) + pushforwardcurl(vals, p) +end + +function (ϕ::GWPCurlRefSpace{T,Deg})(dom::CompScienceMeshes.ReferenceSimplex{Dim}, + u) where {T,Deg,Dim} + + d = Deg + u, v = u + w = 1-u-v + + s = range(zero(T), one(T), length=d+3) + # rwg1 = point(T, u-1, v) + # rwg2 = point(T, u, v-1) + # rwg3 = point(T, u, v) + nd1 = point(T, -v, u-1) + nd2 = point(T, -v+1, u) + nd3 = point(T, -v, u) + + P = SVector{2,T} + vals = P[] + dffu = T[] + dffv = T[] + + i = 0 + for j in 1:d+1 + k = (d+2)-i-j + Rᵢ = _sylpoly_shift(s, i+1, u) + Rⱼ = _sylpoly(s, j+1, v) + Rₖ = _sylpoly(s, k+1, w) + push!(vals, Rᵢ*Rⱼ*Rₖ*nd1) + + # du = _sylpoly_shift_diff(s, i+1, u)*Rⱼ*Rₖ + end + + for i in 1:d+1 + j = 0 + k = (d+2)-i-j + Rᵢ = _sylpoly(s, i+1, u) + Rⱼ = _sylpoly_shift(s, j+1, v) + Rₖ = _sylpoly(s, k+1, w) + push!(vals, Rᵢ*Rⱼ*Rₖ*nd2) + end + + for i in 1:d+1 + j = (d+2)-i + k = 0 + Rᵢ = _sylpoly(s, i+1, u) + Rⱼ = _sylpoly(s, j+1, v) + Rₖ = _sylpoly_shift(s, k+1, w) + push!(vals, Rᵢ*Rⱼ*Rₖ*nd3) + end + + for i in 1:d+1 + for j in 1:d+1 + k = (d+2)-i-j + k <= 0 && continue + Rsᵢ = _sylpoly_shift(s, i+1, u) + Rsⱼ = _sylpoly_shift(s, j+1, v) + Rᵢ = _sylpoly(s, i+1, u) + Rⱼ = _sylpoly(s, j+1, v) + Rₖ = _sylpoly(s, k+1, w) + push!(vals, Rsᵢ*Rⱼ*Rₖ*nd1) + push!(vals, Rᵢ*Rsⱼ*Rₖ*nd2) + end end + + NF = length(vals) + SVector{NF}([(value=f, curl=zero(T)) for f in vals]) +end + + diff --git a/src/bases/local/laglocal.jl b/src/bases/local/laglocal.jl index 63808052..cd90f2c3 100644 --- a/src/bases/local/laglocal.jl +++ b/src/bases/local/laglocal.jl @@ -1,13 +1,13 @@ # T: coeff type # Degree: degree # Dim1: dimension of the support + 1 -struct LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T,NF} end +struct LagrangeRefSpace{T,Degree,Dim1,NF} <: RefSpace{T} end -numfunctions(s::LagrangeRefSpace{T,D,2}) where {T,D} = D+1 -numfunctions(s::LagrangeRefSpace{T,0,3}) where {T} = 1 -numfunctions(s::LagrangeRefSpace{T,1,3}) where {T} = 3 -numfunctions(s::LagrangeRefSpace{T,2,3}) where {T} = 6 -numfunctions(s::LagrangeRefSpace{T,Dg,D1,NF}) where {T,Dg,D1,NF} = NF +numfunctions(s::LagrangeRefSpace{T,D,2}, ch::CompScienceMeshes.ReferenceSimplex{1}) where {T,D} = D+1 +numfunctions(s::LagrangeRefSpace{T,0,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 1 +numfunctions(s::LagrangeRefSpace{T,1,3}, ch::CompScienceMeshes.ReferenceSimplex{2}) where {T} = 3 +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}} @@ -110,7 +110,8 @@ function strace(x::LagrangeRefSpace, cell, localid, face) vals1 = x(P1) vals2 = x(P2) - for j in 1:numfunctions(x) + num_shapes = numfunctions(x, domain(cell)) + for j in 1:num_shapes Q[1,j] = vals1[j].value Q[2,j] = vals2[j].value end @@ -137,12 +138,13 @@ end function restrict(refs::LagrangeRefSpace{T,0}, dom1, dom2) where T - Q = Matrix{T}(I, numfunctions(refs), numfunctions(refs)) + n = numfunctions(refs, domain(dom1)) + Q = Matrix{T}(I, n, n) end function restrict(f::LagrangeRefSpace{T,1}, dom1, dom2) where T - D = numfunctions(f) + D = numfunctions(f, domain(dom1)) Q = zeros(T, D, D) # for each point of the new domain diff --git a/src/bases/local/ncrossbdmlocal.jl b/src/bases/local/ncrossbdmlocal.jl index b7dc7504..5f7fc211 100644 --- a/src/bases/local/ncrossbdmlocal.jl +++ b/src/bases/local/ncrossbdmlocal.jl @@ -1,4 +1,4 @@ -struct NCrossBDMRefSpace{T} <: RefSpace{T,6} end +struct NCrossBDMRefSpace{T} <: RefSpace{T} end function (f::NCrossBDMRefSpace{T})(p) where T diff --git a/src/bases/local/nd2local.jl b/src/bases/local/nd2local.jl index 2a250df0..31dd94c0 100644 --- a/src/bases/local/nd2local.jl +++ b/src/bases/local/nd2local.jl @@ -1,5 +1,5 @@ -mutable struct ND2RefSpace{T} <: RefSpace{T,8} end +mutable struct ND2RefSpace{T} <: RefSpace{T} end function (ϕ::ND2RefSpace)(nbd) @@ -26,3 +26,5 @@ function (ϕ::ND2RefSpace)(nbd) )) end + +numfunctions(x::ND2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8 \ No newline at end of file diff --git a/src/bases/local/ndlcclocal.jl b/src/bases/local/ndlcclocal.jl index ca18947b..37edb79c 100644 --- a/src/bases/local/ndlcclocal.jl +++ b/src/bases/local/ndlcclocal.jl @@ -1,4 +1,4 @@ -struct NDLCCRefSpace{T} <: RefSpace{T,6} end +struct NDLCCRefSpace{T} <: RefSpace{T} end function valuetype(ref::NDLCCRefSpace{T}, charttype::Type) where {T} SVector{universedimension(charttype),T} @@ -58,6 +58,8 @@ function (ϕ::NDLCCRefSpace)(ndlc) ## )) end +numfunctions(x::NDLCCRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 6 + #check orientation function curl(ref::NDLCCRefSpace, sh, el) a = [4,2,3,4,1,2]##[2,1,1,1,2,4]#[2,1,4,4,3,2]#[4,2,3,4,1,2] @@ -71,7 +73,7 @@ end function restrict(ϕ::NDLCCRefSpace{T}, dom1, dom2) where {T} # dom2 is the smaller of the domains - K = numfunctions(ϕ) + K = numfunctions(ϕ, domain(dom1)) D = dimension(dom1) @assert K == 6 diff --git a/src/bases/local/ndlcdlocal.jl b/src/bases/local/ndlcdlocal.jl index ba4cf8ae..1ecc2ad6 100644 --- a/src/bases/local/ndlcdlocal.jl +++ b/src/bases/local/ndlcdlocal.jl @@ -1,4 +1,4 @@ -struct NDLCDRefSpace{T} <: RefSpace{T,4} end +struct NDLCDRefSpace{T} <: RefSpace{T} end function valuetype(ref::NDLCDRefSpace{T}, charttype::Type) where {T} SVector{universedimension(charttype),T} @@ -30,6 +30,8 @@ function (ϕ::NDLCDRefSpace)(ndlc) )) end +numfunctions(x::NDLCDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{3}) = 4 + function ntrace(x::NDLCDRefSpace, el, q, fc) t = zeros(scalartype(x),1,4) t[q] = 1 / volume(fc) @@ -82,7 +84,7 @@ divergence(ref::NDLCDRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el) function restrict(ϕ::NDLCDRefSpace{T}, dom1, dom2) where {T} # dom2 is the smaller of the domains - K = numfunctions(ϕ) + K = numfunctions(ϕ, domain(dom1)) D = dimension(dom1) @assert K == 4 diff --git a/src/bases/local/ndlocal.jl b/src/bases/local/ndlocal.jl index 9f0c196f..a4b3bd17 100644 --- a/src/bases/local/ndlocal.jl +++ b/src/bases/local/ndlocal.jl @@ -6,7 +6,7 @@ This is not the edge starting at vertex `r`. The downside of this local numberin scheme is that it cannot be extended to cells that are not simplices because there is no well defined concept of adjacent-ness. """ -mutable struct NDRefSpace{T} <: RefSpace{T,3} end +mutable struct NDRefSpace{T} <: RefSpace{T} end function (ϕ::NDRefSpace)(nbd) @@ -27,10 +27,11 @@ function (ϕ::NDRefSpace)(nbd) end +numfunctions(x::NDRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3 function restrict(ϕ::NDRefSpace{T}, dom1, dom2) where T - K = numfunctions(ϕ) + K = numfunctions(ϕ, domain(dom1)) D = dimension(dom1) @assert K == 3 diff --git a/src/bases/local/rt2local.jl b/src/bases/local/rt2local.jl index aeaa2c57..bbecede9 100644 --- a/src/bases/local/rt2local.jl +++ b/src/bases/local/rt2local.jl @@ -1,4 +1,4 @@ -struct RT2RefSpace{T} <: RefSpace{T,8} end +struct RT2RefSpace{T} <: RefSpace{T} end function (f::RT2RefSpace)(p) @@ -33,6 +33,7 @@ function (f::RT2RefSpace)(p) ) end +numfunctions(x::RT2RefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 8 function interpolate(fields, interpolant::BEAST.RT2RefSpace, chart) diff --git a/src/bases/local/rtlocal.jl b/src/bases/local/rtlocal.jl index 61e19eac..7c16d6f8 100644 --- a/src/bases/local/rtlocal.jl +++ b/src/bases/local/rtlocal.jl @@ -1,4 +1,4 @@ -struct RTRefSpace{T} <: DivRefSpace{T,3} end +struct RTRefSpace{T} <: DivRefSpace{T} end # valuetype(ref::RTRefSpace{T}, charttype) where {T} = SVector{3,Tuple{SVector{universedimension(charttype),T},T}} function valuetype(ref::RTRefSpace{T}, charttype::Type) where {T} @@ -26,6 +26,8 @@ function (ϕ::RTRefSpace)(mp) )) end +numfunctions(x::RTRefSpace, dom::CompScienceMeshes.ReferenceSimplex{2}) = 3 + divergence(ref::RTRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el)) """ diff --git a/src/bases/local/rtqlocal.jl b/src/bases/local/rtqlocal.jl index bc196d95..f7e38347 100644 --- a/src/bases/local/rtqlocal.jl +++ b/src/bases/local/rtqlocal.jl @@ -1,4 +1,4 @@ -struct RTQuadRefSpace{T} <: DivRefSpace{T,4} end +struct RTQuadRefSpace{T} <: DivRefSpace{T} end function (ϕ::RTQuadRefSpace{T})(p) where {T} diff --git a/src/bases/localbasis.jl b/src/bases/localbasis.jl new file mode 100644 index 00000000..33110983 --- /dev/null +++ b/src/bases/localbasis.jl @@ -0,0 +1,29 @@ +abstract type RefSpace{T} end +abstract type DivRefSpace{T} <: RefSpace{T} end + +function pushforwardcurl(vals, nbhd) + tu = tangents(nbhd, 1) + tv = tangents(nbhd, 2) + j = jacobian(nbhd) + n = normal(nbhd) + map(vals) do v + f = v.value + σ = v.curl + gu = +f[2] + gv = -f[1] + (value=cross(n, (gu*tu + gv*tv)/j), curl=σ/j) + end +end + + +function pushforwarddiv(vals, nbhd) + tu = tangents(nbhd, 1) + tv = tangents(nbhd, 2) + j = jacobian(nbhd) + n = normal(nbhd) + map(vals) do v + f = v.value + σ = v.curl + (value=(f[1]*tu + f[2]*tv)/j, curl=σ/j) + end +end \ No newline at end of file diff --git a/src/bases/subdbasis.jl b/src/bases/subdbasis.jl index bda73d31..9556cd47 100644 --- a/src/bases/subdbasis.jl +++ b/src/bases/subdbasis.jl @@ -14,7 +14,7 @@ function subset(s::S,I) where {S<:subdBasis} end -mutable struct subReferenceSpace{T,D} <: RefSpace{T,D} end +mutable struct subReferenceSpace{T,D} <: RefSpace{T} end refspace(s::subdBasis{T,M,P}) where {T,M,P} = subReferenceSpace{T,12}() # 12 only for regular case diff --git a/src/bases/timebasis.jl b/src/bases/timebasis.jl index dedd754e..4f7fd1ed 100644 --- a/src/bases/timebasis.jl +++ b/src/bases/timebasis.jl @@ -2,7 +2,7 @@ using Compat -mutable struct MonomialBasis{T,Degree,NF} <: RefSpace{T,NF} end +mutable struct MonomialBasis{T,Degree,NF} <: RefSpace{T} end valuetype(::MonomialBasis{T}) where {T} = T degree(::MonomialBasis{T,D}) where {T,D} = D @@ -66,7 +66,7 @@ mutable struct TimeBasisDelta{T} <: AbstractTimeBasisFunction numfunctions::Int end -mutable struct DiracBoundary{T} <: RefSpace{T,1} end +mutable struct DiracBoundary{T} <: RefSpace{T} end numfunctions(x::DiracBoundary) = 1 scalartype(x::TimeBasisDelta{T}) where {T} = T diff --git a/src/excitation.jl b/src/excitation.jl index 165894a7..f194a86c 100644 --- a/src/excitation.jl +++ b/src/excitation.jl @@ -42,13 +42,17 @@ function assemble!(field::Functional, tfs::Space, store; trefs = refspace(tfs) qd = quaddata(field, trefs, tels, quadstrat) + tgeo = geometry(tfs) + tdom = domain(chart(tgeo, first(tgeo))) + num_trefs = numfunctions(trefs, tdom) + for (t, tcell) in enumerate(tels) # compute the testing with the reference elements qr = quadrule(field, trefs, t, tcell, qd, quadstrat) blocal = celltestvalues(trefs, tcell, field, qr) - for i in 1 : numfunctions(trefs) + for i in 1 : num_trefs for (m,a) in tad[t,i] store(a*blocal[i], m) end @@ -82,9 +86,9 @@ function assemble!(field::Functional, tfs::subdBasis, store; end -function celltestvalues(tshs::RefSpace{T,NF}, tcell, field, qr) where {T,NF} +function celltestvalues(tshs::RefSpace{T}, tcell, field, qr) where {T,NF} - num_tshs = numfunctions(tshs) + num_tshs = numfunctions(tshs, domain(tcell)) interactions = zeros(Complex{T}, num_tshs) num_oqp = length(qr) diff --git a/src/helmholtz3d/timedomain/tdhh3dops.jl b/src/helmholtz3d/timedomain/tdhh3dops.jl index afec459a..26911f9b 100644 --- a/src/helmholtz3d/timedomain/tdhh3dops.jl +++ b/src/helmholtz3d/timedomain/tdhh3dops.jl @@ -92,8 +92,8 @@ function innerintegrals!(zlocal, operator::HH3DSingleLayerTDBIO, a = dx / (4*pi) D = operator.num_diffs @assert D == 0 - @assert numfunctions(test_local_space) == 1 - @assert numfunctions(trial_local_space) == 1 + @assert numfunctions(test_local_space, domain(test_element)) == 1 + @assert numfunctions(trial_local_space, domain(trial_element)) == 1 @inline function tmRoR_sl(d, iG) sgn = isodd(d) ? -1 : 1 @@ -139,8 +139,8 @@ function innerintegrals!(zlocal, operator::HH3DHyperSingularTDBIO, trial_element[3], x, r, R, Val{N-1},quad_rule.workspace) - @assert numfunctions(test_local_space) <= 3 - @assert numfunctions(trial_local_space) == 3 + @assert numfunctions(test_local_space, domain(test_element)) <= 3 + @assert numfunctions(trial_local_space, domain(trial_element)) == 3 @inline function tmRoR(d, iG) r = (isodd(d) ? -1 : 1) * iG[d+2] @@ -156,9 +156,9 @@ function innerintegrals!(zlocal, operator::HH3DHyperSingularTDBIO, # weakly singular term α = dx / (4π) * operator.weight_of_weakly_singular_term Ds = operator.num_diffs_on_weakly_singular_term - for i in 1 : numfunctions(test_local_space) + for i in 1 : numfunctions(test_local_space, domain(test_element)) g, curlg = test_values[i] - for j in 1 : numfunctions(trial_local_space) + for j in 1 : numfunctions(trial_local_space, domain(trial_element)) b = trial_element[j] opp_edge = trial_element[mod1(j+2,3)] - trial_element[mod1(j+1,3)] h = norm(opp_edge)/2/volume(trial_element) @@ -175,9 +175,9 @@ function innerintegrals!(zlocal, operator::HH3DHyperSingularTDBIO, # Hyper-singular term β = dx / (4π) * operator.weight_of_hyper_singular_term Dh = operator.num_diffs_on_hyper_singular_term - for i in 1 : numfunctions(test_local_space) + for i in 1 : numfunctions(test_local_space, domain(test_element)) g, curlg = test_values[i] - for j in 1 : numfunctions(trial_local_space) + for j in 1 : numfunctions(trial_local_space, domain(trial_element)) _, curlf = trial_values[j] for k in 1 : numfunctions(time_local_space) d = k - 1 @@ -196,8 +196,8 @@ function innerintegrals!(zlocal, operator::HH3DDoubleLayerTDBIO, test_element, trial_element, time_element, quad_rule, quad_weight) - @assert numfunctions(test_local_space) <= 3 - @assert numfunctions(trial_local_space) == 1 + @assert numfunctions(test_local_space, domain(test_element)) <= 3 + @assert numfunctions(trial_local_space, domain(trial_element)) == 1 dx = quad_weight x = cartesian(test_point) @@ -231,9 +231,9 @@ function innerintegrals!(zlocal, operator::HH3DDoubleLayerTDBIO, # weakly singular term α = dx / (4π) * operator.weight D = operator.num_diffs - for i in 1 : numfunctions(test_local_space) + for i in 1 : numfunctions(test_local_space, domain(test_element)) g, curlg = test_values[i] - for j in 1 : numfunctions(trial_local_space) + for j in 1 : numfunctions(trial_local_space, domain(trial_element)) f, curlf = trial_values[j] for k in 1 : numfunctions(time_local_space) d = k-1 diff --git a/src/helmholtz3d/wiltonints.jl b/src/helmholtz3d/wiltonints.jl index 0921f00a..28772ce1 100644 --- a/src/helmholtz3d/wiltonints.jl +++ b/src/helmholtz3d/wiltonints.jl @@ -52,6 +52,8 @@ function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood, s1, s2, s3 = trial_element.vertices t1, t2, t3 = test_elements.vertices + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) x = cartesian(test_neighborhood) n = normalize((s1-s3)×(s2-s3)) @@ -60,10 +62,10 @@ function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood, scal, vec = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1}) ∫G = (scal[2] + 0.5*γ^2*scal[4]) / (4π) Atot = 1/2*norm(cross(t3-t1,t3-t2)) - for i in 1:numfunctions(test_refspace) + for i in 1:num_tshapes Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x)) g = Ai/Atot - for j in 1:numfunctions(trial_refspace) + for j in 1:num_bshapes zlocal[i,j] += α * ∫G * g * dx end end @@ -81,6 +83,9 @@ function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood, γ = gamma(op) α = op.alpha + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) + s1, s2, s3 = trial_element.vertices x = cartesian(test_neighborhood) @@ -90,7 +95,7 @@ function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood, ∫Rⁿ, ∫RⁿN = WiltonInts84.higherorder(s1,s2,s3,x,3) ∫G = (∫RⁿN[2] + 0.5*γ^2*∫RⁿN[3]) / (4π) - for i in 1:numfunctions(trial_refspace) + for i in 1:num_bshapes zlocal[1,i] += α * ∫G[i] * dx end @@ -110,6 +115,9 @@ function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood, s1, s2, s3 = trial_element.vertices t1, t2, t3 = test_elements.vertices + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) + x = cartesian(test_neighborhood) n = normalize((s1-s3)×(s2-s3)) ρ = x - dot(x - s1, n) * n @@ -118,10 +126,10 @@ function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood, ∫G = (∫RⁿN[2] + 0.5*γ^2*∫RⁿN[3]) / (4π) Atot = 1/2*norm(cross(t3-t1,t3-t2)) - for i in 1:numfunctions(test_refspace) + for i in 1:num_tshapes Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x)) g = Ai/Atot - for j in 1:numfunctions(trial_refspace) + for j in 1:num_bshapes zlocal[i,j] += α * ∫G[j] * g * dx end end @@ -166,12 +174,15 @@ function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood, n = normalize((t1-t3)×(t2-t3)) ρ = x - dot(x - s1, n) * n + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) + scal, vec, grad = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1}) ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π) ∫n∇G = dot(n,∫∇G) - for i in 1:numfunctions(test_refspace) - for j in 1:numfunctions(trial_refspace) + for i in 1:num_tshapes + for j in 1:num_bshapes zlocal[i,j] += α * ∫n∇G * dx end end @@ -194,15 +205,18 @@ function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood, n = normalize((t1-t3)×(t2-t3)) ρ = x - dot(x - s1, n) * n + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) + scal, vec, grad = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1}) ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π) ∫n∇G = dot(n,∫∇G) Atot = 1/2*norm(cross(t3-t1,t3-t2)) - for i in 1:numfunctions(test_refspace) + for i in 1:num_tshapes Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x)) g = Ai/Atot - for j in 1:numfunctions(trial_refspace) + for j in 1:num_bshapes zlocal[i,j] += α * ∫n∇G * g * dx end end @@ -225,11 +239,14 @@ function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood, n = normalize((t1-t3)×(t2-t3)) ρ = x - dot(x - s1, n) * n + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) + _, _, _, grad = WiltonInts84.higherorder(s1,s2,s3,x,3) ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π) - for i in 1:numfunctions(test_refspace) - for j in 1:numfunctions(trial_refspace) + for i in 1:num_tshapes + for j in 1:num_bshapes ∫n∇G = dot(n,∫∇G[j]) zlocal[i,j] += α * ∫n∇G * dx end @@ -253,15 +270,18 @@ function innerintegrals!(op::HH3DDoubleLayerTransposedSng, test_neighborhood, n = normalize((t1-t3)×(t2-t3)) ρ = x - dot(x - s1, n) * n + num_tshapes = numfunctions(test_refspace, domain(test_elements)) + num_bshapes = numfunctions(trial_refspace, domain(trial_element)) + _, _, _, grad = WiltonInts84.higherorder(s1,s2,s3,x,3) ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π) Atot = 1/2*norm(cross(t3-t1,t3-t2)) - for i in 1:numfunctions(test_refspace) + for i in 1:num_tshapes Ai = 1/2*norm(cross(test_elements.vertices[mod1(i-1,3)]-x,test_elements.vertices[mod1(i+1,3)]-x)) g = Ai/Atot - for j in 1:numfunctions(trial_refspace) + for j in 1:num_bshapes ∫n∇G = dot(n,∫∇G[j]) zlocal[i,j] += α * ∫n∇G * g * dx end @@ -303,6 +323,9 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, s1, s2, s3 = s.vertices t1, t2, t3 = t.vertices + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) + x = cartesian(p) n = normalize((s1-s3)×(s2-s3)) ρ = x - dot(x - s1, n) * n @@ -312,10 +335,10 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π) Atot = 1/2*norm(cross(t3-t1,t3-t2)) - for i in 1:numfunctions(g) + for i in 1:num_tshapes Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x)) g = Ai/Atot - for j in 1:numfunctions(f) + for j in 1:num_bshapes z[i,j] += α * dot(n,-∫∇G[j]) * g * dx end end @@ -332,6 +355,9 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, γ = gamma(op) α = op.alpha + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) + s1, s2, s3 = s.vertices t1, t2, t3 = t.vertices @@ -343,8 +369,8 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π) - for i in 1:numfunctions(g) - for j in 1:numfunctions(f) + for i in 1:num_tshapes + for j in 1:num_bshapes z[i,j] += α * dot(n,-∫∇G) * dx #why the minus? end end @@ -363,6 +389,9 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, s1, s2, s3 = s.vertices t1, t2, t3 = t.vertices + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) + x = cartesian(p) n = normalize((s1-s3)×(s2-s3)) ρ = x - dot(x - s1, n) * n @@ -372,13 +401,13 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, ∫∇G = -(grad[1]+0.5*γ^2*grad[3])/(4π) Atot = 1/2*norm(cross(t3-t1,t3-t2)) -for i in 1:numfunctions(g) - Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x)) - g = Ai/Atot - for j in 1:numfunctions(f) - z[i,j] += α * dot(n,-∫∇G) * g * dx + for i in 1:num_tshapes + Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x)) + g = Ai/Atot + for j in 1:num_bshapes + z[i,j] += α * dot(n,-∫∇G) * g * dx + end end -end return nothing end @@ -402,12 +431,14 @@ function innerintegrals!(op::HH3DDoubleLayerSng, p, ∫∇G = -(grad[1] + 0.5*γ^2*grad[2]) / (4π) + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) -for i in 1:numfunctions(g) - for j in 1:numfunctions(f) - z[i,j] += α * dot(n,-∫∇G[j]) * dx + for i in 1:num_tshapes + for j in 1:num_bshapes + z[i,j] += α * dot(n,-∫∇G[j]) * dx + end end -end return nothing end @@ -451,15 +482,18 @@ function innerintegrals!(op::HH3DHyperSingularSng, p, greenconst = (∫Rⁿ[2] + 0.5*γ^2*∫Rⁿ[3]) / (4π) greenlinear = (∫RⁿN[2] + 0.5*γ^2*∫RⁿN[3] ) / (4π) + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) + jt = volume(t) * factorial(dimension(t)) js = volume(s) * factorial(dimension(s)) curlt = [(t3-t2)/jt,(t1-t3)/jt,(t2-t1)/jt] curls = [(s3-s2)/js,(s1-s3)/js,(s2-s1)/js] Atot = 1/2*norm(cross(t3-t1,t3-t2)) - for i in 1:numfunctions(g) + for i in 1:num_tshapes Ai = 1/2*norm(cross(t.vertices[mod1(i-1,3)]-x,t.vertices[mod1(i+1,3)]-x)) g = Ai/Atot - for j in 1:numfunctions(f) + for j in 1:num_bshapes z[i,j] += β * dot(curlt[i],curls[j])*greenconst*dx + α*dot(nx,ny) * greenlinear[j]*g*dx end end diff --git a/src/integralop.jl b/src/integralop.jl index 2989812a..32621123 100644 --- a/src/integralop.jl +++ b/src/integralop.jl @@ -75,12 +75,15 @@ function assemblechunk!(biop::IntegralOperator, tfs::Space, bfs::Space, store; test_elements, tad, tcells = tr bsis_elements, bad, bcells = br - tshapes = refspace(tfs); num_tshapes = numfunctions(tshapes) - bshapes = refspace(bfs); num_bshapes = numfunctions(bshapes) - tgeo = geometry(tfs) bgeo = geometry(bfs) + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + tshapes = refspace(tfs); num_tshapes = numfunctions(tshapes, tdom) + bshapes = refspace(bfs); num_bshapes = numfunctions(bshapes, bdom) + qs = if CompScienceMeshes.refines(tgeo, bgeo) TestRefinesTrialQStrat(quadstrat) elseif CompScienceMeshes.refines(bgeo, tgeo) @@ -328,8 +331,14 @@ function assembleblock_primer(biop, tfs, bfs; test_elements, tad = assemblydata(tfs; onlyactives=false) bsis_elements, bad = assemblydata(bfs; onlyactives=false) - tshapes = refspace(tfs); num_tshapes = numfunctions(tshapes) - bshapes = refspace(bfs); num_bshapes = numfunctions(bshapes) + tgeo = geometry(tfs) + bgeo = geometry(bfs) + + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + tshapes = refspace(tfs); num_tshapes = numfunctions(tshapes, tdom) + bshapes = refspace(bfs); num_bshapes = numfunctions(bshapes, bdom) qd = quaddata(biop, tshapes, bshapes, test_elements, bsis_elements, quadstrat) @@ -551,14 +560,20 @@ end end end end end end end function assemblerow!(biop::IntegralOperator, test_functions::Space, trial_functions::Space, store; quadstrat=defaultquadstrat(biop, test_functions, trial_functions)) - test_elements = elements(geometry(test_functions)) + tgeo = geometry(test_functions) + bgeo = geometry(trial_functions) + + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + test_elements = elements(tgeo) trial_elements, trial_assembly_data = assemblydata(trial_functions) test_shapes = refspace(test_functions) trial_shapes = refspace(trial_functions) - num_test_shapes = numfunctions(test_shapes) - num_trial_shapes = numfunctions(trial_shapes) + num_test_shapes = numfunctions(test_shapes, tdom) + num_trial_shapes = numfunctions(trial_shapes, bdom) quadrature_data = quaddata(biop, test_shapes, trial_shapes, test_elements, trial_elements, quadstrat) diff --git a/src/localop.jl b/src/localop.jl index 9ce1410b..a169ae7d 100644 --- a/src/localop.jl +++ b/src/localop.jl @@ -99,12 +99,21 @@ function assemble_local_matched!(biop::LocalOperator, tfs::Space, bfs::Space, st trefs = refspace(tfs) brefs = refspace(bfs) + tgeo = geometry(tfs) + bgeo = geometry(bfs) + + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + num_trefs = numfunctions(trefs, tdom) + num_brefs = numfunctions(brefs, bdom) + qd = quaddata(biop, trefs, brefs, tels, bels, quadstrat) verbose = length(tels) > 10_000 verbose && print("dots out of 20: ") todo, done, pctg = length(tels), 0, 0 - locmat = zeros(scalartype(biop, trefs, brefs), numfunctions(trefs), numfunctions(brefs)) + locmat = zeros(scalartype(biop, trefs, brefs), num_trefs, num_brefs) for (p,cell) in enumerate(tels) P = ta2g[p] q = bg2a[P] @@ -138,6 +147,15 @@ function assemble_local_refines!(biop::LocalOperator, tfs::Space, bfs::Space, st trefs = refspace(tfs) brefs = refspace(bfs) + tgeo = geometry(tfs) + bgeo = geometry(bfs) + + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + num_trefs = numfunctions(trefs, tdom) + num_brefs = numfunctions(brefs, bdom) + tels, tad, ta2g = assemblydata(tfs) bels, bad, ba2g = assemblydata(bfs) @@ -167,8 +185,8 @@ function assemble_local_refines!(biop::LocalOperator, tfs::Space, bfs::Space, st zlocal = cellinteractions(biop, trefs, brefs, cell, qr) zlocal = Q * zlocal * P' - for i in 1 : numfunctions(trefs) - for j in 1 : numfunctions(brefs) + for i in 1 : num_trefs + for j in 1 : num_brefs for (m,a) in tad[p,i] for (n,b) in bad[q,j] store(a * zlocal[i,j] * b, m, n) @@ -256,6 +274,15 @@ function assemble_local_mixed!(biop::LocalOperator, tfs::Space{T}, bfs::Space{T} trefs = refspace(tfs) brefs = refspace(bfs) + tgeo = geometry(tfs) + bgeo = geometry(bfs) + + tdom = domain(chart(tgeo, first(tgeo))) + bdom = domain(chart(bgeo, first(bgeo))) + + num_trefs = numfunctions(trefs, tdom) + num_brefs = numfunctions(brefs, bdom) + tels, tad = assemblydata(tfs) bels, bad = assemblydata(bfs) @@ -288,8 +315,8 @@ function assemble_local_mixed!(biop::LocalOperator, tfs::Space{T}, bfs::Space{T} zlocal = cellinteractions(biop, trefs, brefs, cell, qr) zlocal = Q * zlocal * P' - for i in 1 : numfunctions(trefs) - for j in 1 : numfunctions(brefs) + for i in 1 : num_trefs + for j in 1 : num_brefs for (m,a) in tad[p,i] for (n,b) in bad[q,j] store(a * zlocal[i,j] * b, m, n) diff --git a/src/maxwell/timedomain/mwtdops.jl b/src/maxwell/timedomain/mwtdops.jl index 962e577f..86f69d2b 100644 --- a/src/maxwell/timedomain/mwtdops.jl +++ b/src/maxwell/timedomain/mwtdops.jl @@ -266,10 +266,13 @@ function innerintegrals!(zl, op::MWSingleLayerTDIO, sol5 = sol4*sol solpowers = (one(sol), sol, sol2, sol3, sol4, sol5) - for i in 1 : numfunctions(U) + udim = numfunctions(U, domain(τ)) + vdim = numfunctions(V, domain(σ)) + + for i in 1 : udim a = τ[i] g = (x-a) - for j in 1 : numfunctions(V) + for j in 1 : vdim b = σ[j]; bξ = ξ-b for k in 1 : numfunctions(W) d = k-1 # ranges from 0 to numfunctions(W)-1 @@ -338,11 +341,14 @@ function innerintegrals!(z, op::MWDoubleLayerTDIO, Ux = U(p) Vx = αf * @SVector[(x-σ[1]), (x-σ[2]), (x-σ[3])] - for i in 1 : numfunctions(U) + udim = numfunctions(U, domain(τ)) + vdim = numfunctions(V, domain(σ)) + + for i in 1 : udim # a = τ[i] # g = αg * (x-τ[i]) g = Ux[i].value - for j in 1 : numfunctions(V) + for j in 1 : vdim # b = σ[j] # f = αf * (x-σ[j]) # f = Vx[j].value diff --git a/src/maxwell/wiltonints.jl b/src/maxwell/wiltonints.jl index 4ffe2d80..9b715446 100644 --- a/src/maxwell/wiltonints.jl +++ b/src/maxwell/wiltonints.jl @@ -29,13 +29,16 @@ function innerintegrals!(op::MWSingleLayer3DSng, p, g, f, t, s, z, c₁ = op.α c₂ = op.β + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) + α = 1 / volume(t) / volume(s) / 4 - for i in 1 : numfunctions(g) + for i in 1 : num_bshapes a = t[i] g = x - a dg = 2 - for j in 1 : numfunctions(f) + for j in 1 : num_tshapes b = s[j] ∫Gf = SVector(∫Gy[1]-∫G*b[1], ∫Gy[2]-∫G*b[2], ∫Gy[3]-∫G*b[3]) @@ -60,14 +63,17 @@ function innerintegrals!(op::MWDoubleLayer3DSng, p, g, f, t, s, z, strat::Wilton scal, vec, grad = WiltonInts84.wiltonints(s[1], s[2], s[3], x, Val{1}) + num_tshapes = numfunctions(g, domain(t)) + num_bshapes = numfunctions(f, domain(s)) + # \int \nabla G_s with G_s = \nabla (1/R + 0.5*γ^2*R) / (4\pi) ∫∇G = T.((-grad[1] - 0.5*γ^2*grad[3]) / (4π)) α = 1 / volume(t) / volume(s) / 4 - for i in 1 : numfunctions(g) + for i in 1 : num_tshapes a = t[i] g = (x - a) - for j in 1 : numfunctions(f) + for j in 1 : num_bshapes b = s[j] z[i,j] += ( α * ( (x-b) × g ) ⋅ ∫∇G ) * dx diff --git a/src/postproc.jl b/src/postproc.jl index 961ad2e9..398f4c8f 100644 --- a/src/postproc.jl +++ b/src/postproc.jl @@ -171,7 +171,9 @@ function potential!(store, op, points, basis; els, ad = assemblydata(basis) rs = refspace(basis) - zlocal = Array{type}(undef,numfunctions(rs)) + geo = geometry(basis) + nf = numfunctions(rs, domain(chart(geo, first(geo)))) + zlocal = Array{type}(undef, nf) qdata = quaddata(op,rs,els,quadstrat) print("dots out of 10: ") diff --git a/src/quadrature/rules/testrefinestrialqrule.jl b/src/quadrature/rules/testrefinestrialqrule.jl index 08bcec5a..e2d86b99 100644 --- a/src/quadrature/rules/testrefinestrialqrule.jl +++ b/src/quadrature/rules/testrefinestrialqrule.jl @@ -13,6 +13,12 @@ function momintegrals!(out, op, test_mesh = geometry(test_functions) trial_mesh = geometry(trial_functions) + tdom = domain(test_chart) + bdom = domain(trial_chart) + + num_tshapes = numfunctions(test_local_space, tdom) + num_bshapes = numfunctions(trial_local_space, bdom) + parent_mesh = CompScienceMeshes.parent(test_mesh) trial_charts = [chart(test_mesh, p) for p in CompScienceMeshes.children(parent_mesh, trial_cell)] @@ -31,8 +37,8 @@ function momintegrals!(out, op, test_functions, nothing, test_chart, trial_functions, nothing, chart, qr) - for j in 1:numfunctions(trial_local_space) - for i in 1:numfunctions(test_local_space) + for j in 1:num_bshapes + for i in 1:num_tshapes for k in 1:size(Q, 2) out[i,j] += zlocal[i,k] * Q[j,k] end end end end end \ No newline at end of file diff --git a/src/quadrature/rules/trialrefinestestqrule.jl b/src/quadrature/rules/trialrefinestestqrule.jl index 3bbb4c7a..f0c74f54 100644 --- a/src/quadrature/rules/trialrefinestestqrule.jl +++ b/src/quadrature/rules/trialrefinestestqrule.jl @@ -13,6 +13,12 @@ function momintegrals!(out, op, test_mesh = geometry(test_functions) trial_mesh = geometry(trial_functions) + tdom = domain(chart(test_mesh, first(test_mesh))) + bdom = domain(chart(trial_mesh, first(trial_mesh))) + + num_tshapes = numfunctions(test_local_space, tdom) + num_bshapes = numfunctions(trial_local_space, bdom) + parent_mesh = CompScienceMeshes.parent(trial_mesh) test_charts = [chart(trial_mesh, p) for p in CompScienceMeshes.children(parent_mesh, test_cell)] @@ -30,8 +36,8 @@ function momintegrals!(out, op, test_functions, nothing, chart, trial_functions, trial_cell, trial_chart, qr) - for j in 1:numfunctions(trial_local_space) - for i in 1:numfunctions(test_local_space) + for j in 1:num_bshapes + for i in 1:num_tshapes for k in 1:size(Q, 2) out[i,j] += Q[i,k] * zlocal[k,j] end end end end end \ No newline at end of file diff --git a/src/quadrature/sauterschwabints.jl b/src/quadrature/sauterschwabints.jl index 423bda9a..1d56cdb3 100644 --- a/src/quadrature/sauterschwabints.jl +++ b/src/quadrature/sauterschwabints.jl @@ -150,10 +150,13 @@ function momintegrals!(op::Operator, vertices(test_chart), vertices(trial_chart), rule) + num_tshapes = numfunctions(test_local_space, domain(test_chart)) + num_bshapes = numfunctions(trial_local_space, domain(trial_chart)) + igd = Integrand(op, test_local_space, trial_local_space, test_chart, trial_chart) igdp = pulledback_integrand(igd, I, test_chart, J, trial_chart) G = SauterSchwabQuadrature.sauterschwab_parameterized(igdp, rule) - out[1:numfunctions(test_local_space),1:numfunctions(trial_local_space)] .+= G + out[1:num_tshapes, 1:num_bshapes] .+= G nothing end diff --git a/src/timedomain/tdexcitation.jl b/src/timedomain/tdexcitation.jl index cf3819af..e6805273 100644 --- a/src/timedomain/tdexcitation.jl +++ b/src/timedomain/tdexcitation.jl @@ -85,9 +85,11 @@ function assemble!(exc::TDFunctional, testST, store; testels, testad = assemblydata(testfns) timeels, timead = assemblydata(timefns) + qd = quaddata(exc, testrefs, timerefs, testels, timeels) - - z = zeros(eltype(exc), numfunctions(testrefs), numfunctions(timerefs)) + + num_testshapes = numfunctions(testrefs, domain(first(testels))) + z = zeros(eltype(exc), num_testshapes, numfunctions(timerefs)) for p in eachindex(testels) τ = testels[p] for r in eachindex(timeels) @@ -97,7 +99,7 @@ function assemble!(exc::TDFunctional, testST, store; qr = quadrule(exc, testrefs, timerefs, p, τ, r, ρ, qd) momintegrals!(z, exc, testrefs, timerefs, τ, ρ, qr) - for i in 1 : numfunctions(testrefs) + for i in 1 : num_testshapes for d in 1 : numfunctions(timerefs) v = z[i,d] @@ -155,13 +157,15 @@ end function timeintegrals!(z, exc::TDFunctional, testrefs, timerefs, testpoint, timeelement, dx, qr, f) + num_tshapes = numfunctions(testrefs, domain(chart(testpoint))) + for p in qr.quad_points t = p.point w = p.weight U = p.value dt = w #* jacobian(t) # * volume(timeelement) - for i in 1 : numfunctions(testrefs) + for i in 1 : num_tshapes for k in 1 : numfunctions(timerefs) z[i,k] += dot(f[i][1]*U[k], exc(testpoint,t)) * dt * dx end @@ -176,12 +180,14 @@ function timeintegrals!(z, exc::TDFunctional, testpoint, timeelement, dx, qr, testvals) + num_tshapes = numfunctions(spacerefs, domain(chart(testpoint))) + # since timeelement uses barycentric coordinates, # the first/left vertex has coords u = 1.0! testtime = neighborhood(timeelement, point(0.0)) @assert cartesian(testtime)[1] ≈ timeelement.vertices[2][1] - for i in 1 : numfunctions(spacerefs) + for i in 1 : num_tshapes z[i,1] += dot(testvals[i][1], exc(testpoint, testtime)) * dx end end diff --git a/src/timedomain/tdintegralop.jl b/src/timedomain/tdintegralop.jl index bd9820f6..73bcb201 100644 --- a/src/timedomain/tdintegralop.jl +++ b/src/timedomain/tdintegralop.jl @@ -251,8 +251,11 @@ function assemble_chunk!(op::RetardedPotential, testST, trialST, store; qd = quaddata(op, U, V, W, testels, trialels, nothing, quadstrat) - udim = numfunctions(U) - vdim = numfunctions(V) + ugeo = geometry(testspace) + vgeo = geometry(trialspace) + + udim = numfunctions(U, domain(chart(ugeo, first(ugeo)))) + vdim = numfunctions(V, domain(chart(vgeo, first(vgeo)))) wdim = numfunctions(W) z = zeros(scalartype(op, testST, trialST), udim, vdim, wdim) diff --git a/src/utils/lagpolys.jl b/src/utils/lagpolys.jl new file mode 100644 index 00000000..65d3c924 --- /dev/null +++ b/src/utils/lagpolys.jl @@ -0,0 +1,45 @@ +function _lagpoly(nodes, i, s, i0=1, i1=length(nodes)) + T = eltype(s) + r = one(T) + si = nodes[i] + for j in i0:i1 + j == i && continue + sj = nodes[j] + r *= (s - sj) / (si - sj) + end + return r +end + +function _lagpoly_diff(nodes, i, s, i0=1, i1=length(nodes)) + r = zero(T) + si = nodes[i] + for p in i0:i1 + p == i && continue + rp = one(T) + for j in i0:i1 + j == i && continue + j == p && continue + sj = nodees[j] + rp *= (s - sj) / (si - sj) + end + rp *= 1 / (si - sp) + r += rp + end + return r +end + +function _sylpoly(nodes, i, s) + _lagpoly(nodes, i, s, 1, i) +end + +function _sylpoly_diff(nodes, i, s) + _lagpoly_diff(nodes, i, s, 1, i) +end + +function _sylpoly_shift(nodes, i, s) + _lagpoly(nodes, i, s, 2, i) +end + +function _sylpoly_shift_diff(nodes, i, s) + _lagpoly_shift(nodes, i, s, 2, i) +end \ No newline at end of file diff --git a/test/test_basis.jl b/test/test_basis.jl index c6f36017..a85c7937 100644 --- a/test/test_basis.jl +++ b/test/test_basis.jl @@ -75,8 +75,9 @@ for T in [Float32, Float64] m = meshrectangle(T(1.0), T(1.0), T(0.5), 3) X = lagrangec0d1(m) x = refspace(X) + dom = domain(chart(m, first(m))) - @test numfunctions(x) == 3 + @test numfunctions(x, dom) == 3 @test numfunctions(X) == 1 @test length(X.fns[1]) == 6 @@ -213,9 +214,10 @@ x = refspace(X) isonjunction = inclosure_gpredicate(j) els, ad = BEAST.assemblydata(X) +num_shapes = numfunctions(x, domain(chart(m, first(m)))) for _p in 1:numcells(m) el = els[_p] - for r in 1:numfunctions(x) + for r in 1:num_shapes vert = el[r] isonjunction(vert) || continue for (i,w) in ad[_p,r] diff --git a/test/test_gwp.jl b/test/test_gwp.jl new file mode 100644 index 00000000..7d1b83b7 --- /dev/null +++ b/test/test_gwp.jl @@ -0,0 +1,12 @@ +@testitem "GWPRefSpace eval" begin + using CompScienceMeshes + + T = Float64 + Dim, Nverts, Degree = 2, 3, 3 + + dom = CompScienceMeshes.ReferenceSimplex{Dim,T,Nverts}() + ϕ = BEAST.GWPCurlRefSpace{T,Degree}() + + u = (one(T)/3, one(T)/3) + vals = ϕ(dom, u) +end \ No newline at end of file diff --git a/test/test_rt.jl b/test/test_rt.jl index 56d9b781..ea46e92c 100644 --- a/test/test_rt.jl +++ b/test/test_rt.jl @@ -13,7 +13,8 @@ tol = eps(T) * 10^3 function shapevals(ϕ, ts) numpoints = length(ts) - numshapes = numfunctions(ϕ) + dom = CompScienceMeshes.domain(CompScienceMeshes.chart(ts[1])) + numshapes = numfunctions(ϕ, dom) y = Array{P}(undef, numshapes, numpoints) for i in 1 : numpoints diff --git a/test/test_ss_nested_meshes.jl b/test/test_ss_nested_meshes.jl index 48bf2c98..b23f725d 100644 --- a/test/test_ss_nested_meshes.jl +++ b/test/test_ss_nested_meshes.jl @@ -54,7 +54,7 @@ end qs_strat = BEAST.DoubleNumWiltonSauterQStrat(1,1,6,7,10,10,10,10) # sauterschwab = BEAST.SauterSchwabQuadrature.CommonFace(BEAST._legendre(10,0.0,1.0)) -out_ss = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(out_ss, 𝒜, Y, p, test_chart, X, 1, trial_chart, @@ -62,7 +62,7 @@ BEAST.momintegrals!(out_ss, 𝒜, wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(out_dw, 𝒜, Y, p, test_chart, X, 1, trial_chart, @@ -77,11 +77,11 @@ test_quadpoints = BEAST.quadpoints(𝒳, [test_chart], (12,))[1,1] trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1] sauterschwab = BEAST.SauterSchwabQuadrature.CommonFace(BEAST._legendre(10,0.0,1.0)) -out_ss1 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss1 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss1,sauterschwab) wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw1 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw1 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw1,wiltonsingext) @test out_ss1 ≈ out_dw1 rtol=3e-6 @@ -94,11 +94,11 @@ test_quadpoints = BEAST.quadpoints(𝒳, [test_chart], (12,))[1,1] trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1] sauterschwab = BEAST.SauterSchwabQuadrature.CommonEdge(BEAST._legendre(10,0.0,1.0)) -out_ss2 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss2 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss2,sauterschwab) wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw2 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw2 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw2,wiltonsingext) @test out_ss2 ≈ out_dw2 rtol=4e-6 @@ -111,11 +111,11 @@ test_quadpoints = BEAST.quadpoints(𝒳, [test_chart], (12,))[1,1] trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1] sauterschwab = BEAST.SauterSchwabQuadrature.CommonEdge(BEAST._legendre(10,0.0,1.0)) -out_ss3 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss3 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss3,sauterschwab) wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw3 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw3 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw3,wiltonsingext) @test out_ss3 ≈ out_dw3 rtol=3e-6 @@ -129,11 +129,11 @@ test_quadpoints = BEAST.quadpoints(𝒳, [test_chart], (12,))[1,1] trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1] sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0)) -out_ss4 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss4 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss4,sauterschwab) wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw4 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw4 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw4,wiltonsingext) @test out_ss4 ≈ out_dw4 rtol=2e-9 @@ -146,11 +146,11 @@ test_quadpoints = BEAST.quadpoints(𝒳, [test_chart], (12,))[1,1] trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1] sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0)) -out_ss5 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss5 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss5,sauterschwab) wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw5 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw5 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw5,wiltonsingext) @test out_ss4 ≈ out_dw4 rtol=3e-8 @@ -163,11 +163,11 @@ test_quadpoints = BEAST.quadpoints(𝒳, [test_chart], (12,))[1,1] trial_quadpoints = BEAST.quadpoints(𝒳, [trial_chart], (13,))[1,1] sauterschwab = BEAST.SauterSchwabQuadrature.CommonVertex(BEAST._legendre(10,0.0,1.0)) -out_ss6 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_ss6 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_ss6,sauterschwab) wiltonsingext = BEAST.WiltonSERule(test_quadpoints, BEAST.DoubleQuadRule(test_quadpoints, trial_quadpoints)) -out_dw6 = zeros(T, numfunctions(𝒳), numfunctions(𝒳)) +out_dw6 = zeros(T, numfunctions(𝒳, domain(test_chart)), numfunctions(𝒳, domain(trial_chart))) BEAST.momintegrals!(𝒜,𝒳,𝒳,test_chart,trial_chart,out_dw6,wiltonsingext) @test out_ss4 ≈ out_dw4 rtol=6e-8 diff --git a/test/test_tdassembly.jl b/test/test_tdassembly.jl index feac5229..00894b9f 100644 --- a/test/test_tdassembly.jl +++ b/test/test_tdassembly.jl @@ -58,14 +58,14 @@ function momintegrals!(z, op::typeof(G), TR = [(-R)^d for d in 0:maxdegree] dxdy = wpv_τ.weight * wpv_σ.weight - for i in 1 : numfunctions(g) - for j in 1 : numfunctions(f) + for i in 1 : numfunctions(g, domain(τ)) + for j in 1 : numfunctions(f, domain(σ)) for k in 1 : numfunctions(T) z[i,j,k] += dxdy * gx[i] * fy[j] * TR[k] / (4π*R) end end end end end end # Compare results for a single monomial -z1 = zeros(numfunctions(x1), numfunctions(x2), numfunctions(q)) +z1 = zeros(numfunctions(x1, domain(τ1)), numfunctions(x2, domain(τ2)), numfunctions(q)) for r in BEAST.rings(τ1, τ2, ΔR) local ι = BEAST.ring(r, ΔR) momintegrals!(z1, G, x1, x2, q, τ1, τ2, ι, DoubleQuadTimeDomainRule()) @@ -73,7 +73,7 @@ end qs = BEAST.defaultquadstrat(G,X1,X2) qd = quaddata(G, x1, x2, q, [τ1], [τ2], nothing, qs) -z2 = zeros(numfunctions(x1), numfunctions(x2), numfunctions(q)) +z2 = zeros(numfunctions(x1, domain(τ1)), numfunctions(x2, domain(τ2)), numfunctions(q)) for r in BEAST.rings(τ1, τ2, ΔR) local ι = BEAST.ring(r, ΔR) quad_rule = quadrule(G, x1, x2, q, 1, τ1, 1, τ2, r, ι, qd, qs) diff --git a/test/test_tdhhdbl.jl b/test/test_tdhhdbl.jl index f56de5e0..864fd295 100644 --- a/test/test_tdhhdbl.jl +++ b/test/test_tdhhdbl.jl @@ -13,11 +13,13 @@ c = 1.0 K = BEAST.HH3DDoubleLayerTDBIO(speed_of_light=c) X = lagrangecxd0(G) -@test numfunctions(refspace(X)) == 1 +@test numfunctions(refspace(X), domain(chart(G, first(G)))) == 1 Y = duallagrangec0d1(G) @test numfunctions(Y) == 1 -@test numfunctions(refspace(Y)) == 3 + +Gref = geometry(Y) +@test numfunctions(refspace(Y), domain(chart(Gref, first(Gref)))) == 3 # X = duallagrangecxd0(G, boundary(G)) # Y = lagrangec0d1(G, dirichlet=false) diff --git a/test/test_wiltonints.jl b/test/test_wiltonints.jl index ed0acfbe..1f1469b7 100644 --- a/test/test_wiltonints.jl +++ b/test/test_wiltonints.jl @@ -333,7 +333,7 @@ x = BEAST.refspace(X) κ = 1.0 op = BEAST.MWSingleLayer3D(κ) -n = BE.numfunctions(x) +n = BE.numfunctions(x, domain(s)) z1 = zeros(ComplexF64, n, n) z2 = zeros(ComplexF64, n, n)