Skip to content

Commit

Permalink
consolidation quadstrats - tests pass
Browse files Browse the repository at this point in the history
  • Loading branch information
krcools committed Jun 28, 2024
1 parent be8a5a4 commit 0caed85
Show file tree
Hide file tree
Showing 21 changed files with 214 additions and 340 deletions.
20 changes: 11 additions & 9 deletions src/BEAST.jl
Original file line number Diff line number Diff line change
Expand Up @@ -185,15 +185,6 @@ include("bases/tensorbasis.jl")
include("operator.jl")

include("quadrature/quadstrats.jl")

include("excitation.jl")
include("localop.jl")
include("multiplicativeop.jl")
include("identityop.jl")
include("integralop.jl")
include("dyadicop.jl")
include("interpolation.jl")

include("quadrature/doublenumqstrat.jl")
include("quadrature/doublenumsauterqstrat.jl")
include("quadrature/doublenumwiltonsauterqstrat.jl")
Expand All @@ -202,12 +193,23 @@ include("quadrature/selfsauterdnumotherwiseqstrat.jl")
include("quadrature/nonconformingintegralopqstrat.jl")
include("quadrature/commonfaceoverlappingedgeqstrat.jl")
include("quadrature/strategies/cfcvsautercewiltonpdnumqstrat.jl")
include("quadrature/strategies/testrefinestrialqstrat.jl")

include("excitation.jl")
include("localop.jl")
include("multiplicativeop.jl")
include("identityop.jl")
include("integralop.jl")
include("dyadicop.jl")
include("interpolation.jl")

include("quadrature/rules/momintegrals.jl")
include("quadrature/doublenumints.jl")
include("quadrature/singularityextractionints.jl")
include("quadrature/sauterschwabints.jl")
include("quadrature/nonconformingoverlapqrule.jl")
include("quadrature/nonconformingtouchqrule.jl")
include("quadrature/rules/testrefinestrialqrule.jl")

include("postproc.jl")
include("postproc/segcurrents.jl")
Expand Down
11 changes: 6 additions & 5 deletions src/decoupled/dpops.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,14 @@ function integrand(op::CurlSingleLayerDP3D, kernel_vals,
return -α * dot(nx × gx, ∇G * fy)
end

function momintegrals!(op::CurlSingleLayerDP3D,
test_local_space::RTRefSpace, trial_local_space::LagrangeRefSpace,
test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)
function momintegrals!(out, op::CurlSingleLayerDP3D,
test_local_space::RTRefSpace, tptr, test_triangular_element,
trial_local_space::LagrangeRefSpace, bptr, trial_triangular_element,
qrule::SauterSchwabStrategy)

I, J, K, L = SauterSchwabQuadrature.reorder(
test_triangular_element.vertices,
trial_triangular_element.vertices, strat)
trial_triangular_element.vertices, qrule)

test_triangular_element = simplex(test_triangular_element.vertices[I]...)
trial_triangular_element = simplex(trial_triangular_element.vertices[J]...)
Expand Down Expand Up @@ -55,7 +56,7 @@ function momintegrals!(op::CurlSingleLayerDP3D,
return @SMatrix[dot(f[i].value × nx, R[j]) for i in 1:3, j in 1:3]
end

Q = sauterschwab_parameterized(igd, strat)
Q = sauterschwab_parameterized(igd, qrule)
for j 1:3
for i 1:3
out[i,j] += Q[K[i],L[j]]
Expand Down
4 changes: 3 additions & 1 deletion src/helmholtz2d/helmholtzop.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,9 @@ function testfunc1()
print("test function!")
end

defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumWiltonSauterQStrat(4,3,4,3,4,4,4,4)
# defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumWiltonSauterQStrat(4,3,4,3,4,4,4,4)
defaultquadstrat(op::HelmholtzOperator2D, tfs, bfs) = DoubleNumQStrat(4,3)


function quaddata(op::HelmholtzOperator2D, g::LagrangeRefSpace, f::LagrangeRefSpace, tels, bels,
qs::DoubleNumWiltonSauterQStrat)
Expand Down
264 changes: 0 additions & 264 deletions src/helmholtz3d/hh3d_sauterschwabqr.jl
Original file line number Diff line number Diff line change
@@ -1,265 +1 @@
# function pulled_back_integrand(op::HH3DSingleLayerFDBIO,
# test_local_space::LagrangeRefSpace,
# trial_local_space::LagrangeRefSpace,
# test_chart, trial_chart)

# (u,v) -> begin

# x = neighborhood(test_chart,u)
# y = neighborhood(trial_chart,v)

# f = test_local_space(x)
# g = trial_local_space(y)

# j = jacobian(x) * jacobian(y)

# α = op.alpha
# γ = gamma(op)
# R = norm(cartesian(x)-cartesian(y))
# G = exp(-γ*R)/(4*π*R)

# αjG = α*G*j

# SMatrix{length(f),length(g)}((f[j].value * αjG * g[i].value for i in 1:length(g) for j in 1:length(f) )...)
# end
# end


# function pulled_back_integrand(op::HH3DHyperSingularFDBIO,
# test_local_space::LagrangeRefSpace,
# trial_local_space::LagrangeRefSpace,
# test_chart, trial_chart)

# (u,v) -> begin

# x = neighborhood(test_chart,u)
# y = neighborhood(trial_chart,v)

# nx = normal(x)
# ny = normal(y)

# f = test_local_space(x)
# g = trial_local_space(y)

# j = jacobian(x) * jacobian(y)

# α = op.alpha
# β = op.beta
# γ = gamma(op)
# R = norm(cartesian(x)-cartesian(y))
# G = exp(-γ*R)/(4*π*R)

# αjG = ny*α*G*j
# βjG = β*G*j

# A = SA[(αjG*g[i].value for i in 1:length(g))...]
# B = SA[(βjG*g[i].curl for i in 1:length(g))...]

# SMatrix{length(f),length(g)}((((dot(nx*f[j].value,A[i])+dot(f[j].curl,B[i])) for i in 1:length(g) for j in 1:length(f))...))
# end
# end

# function pulled_back_integrand(op::HH3DDoubleLayerFDBIO,
# test_local_space::LagrangeRefSpace,
# trial_local_space::LagrangeRefSpace,
# test_chart, trial_chart)

# (u,v) -> begin

# x = neighborhood(test_chart,u)
# y = neighborhood(trial_chart,v)

# ny = normal(y)
# f = test_local_space(x)
# g = trial_local_space(y)

# j = jacobian(x) * jacobian(y)

# α = op.alpha
# γ = gamma(op)

# r = cartesian(x) - cartesian(y)
# R = norm(r)
# G = exp(-γ*R)/(4*π*R)
# inv_R = 1/R
# ∇G = -(γ + inv_R) * G * inv_R * r
# αnyj∇G = dot(ny,-α*∇G*j)

# SMatrix{length(f),length(g)}((f[j].value * αnyj∇G * g[i].value for i in 1:length(g) for j in 1:length(f))...)
# end
# end

# function pulled_back_integrand(op::HH3DDoubleLayerTransposedFDBIO,
# test_local_space::LagrangeRefSpace,
# trial_local_space::LagrangeRefSpace,
# test_chart, trial_chart)

# (u,v) -> begin

# x = neighborhood(test_chart,u)
# y = neighborhood(trial_chart,v)

# nx = normal(x)
# f = test_local_space(x)
# g = trial_local_space(y)

# j = jacobian(x) * jacobian(y)

# α = op.alpha
# γ = gamma(op)

# r = cartesian(x) - cartesian(y)
# R = norm(r)
# G = exp(-γ*R)/(4*π*R)
# inv_R = 1/R
# ∇G = -(γ + inv_R) * G * inv_R * r
# αnxj∇G = dot(nx,α*∇G*j)

# SMatrix{length(f),length(g)}((f[j].value * αnxj∇G * g[i].value for i in 1:length(g) for j in 1:length(f))...)
# end
# end

# function momintegrals!(op::Helmholtz3DOp,
# test_local_space::LagrangeRefSpace{<:Any,0}, trial_local_space::LagrangeRefSpace{<:Any,0},
# test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

# I, J, K, L = SauterSchwabQuadrature.reorder(
# test_triangular_element.vertices,
# trial_triangular_element.vertices, strat)

# test_triangular_element = simplex(
# test_triangular_element.vertices[I[1]],
# test_triangular_element.vertices[I[2]],
# test_triangular_element.vertices[I[3]])

# trial_triangular_element = simplex(
# trial_triangular_element.vertices[J[1]],
# trial_triangular_element.vertices[J[2]],
# trial_triangular_element.vertices[J[3]])

# test_sign = Combinatorics.levicivita(I)
# trial_sign = Combinatorics.levicivita(J)
# σ = momintegrals_sign(op, test_sign, trial_sign)

# igd = pulled_back_integrand(op, test_local_space, trial_local_space,
# test_triangular_element, trial_triangular_element)
# G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)
# out[1,1] += G[1,1] * σ

# nothing
# end


# function momintegrals!(op::Helmholtz3DOp,
# test_local_space::LagrangeRefSpace, trial_local_space::LagrangeRefSpace,
# test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

# I, J, K, L = SauterSchwabQuadrature.reorder(
# test_triangular_element.vertices,
# trial_triangular_element.vertices, strat)

# test_triangular_element = simplex(
# test_triangular_element.vertices[I[1]],
# test_triangular_element.vertices[I[2]],
# test_triangular_element.vertices[I[3]])

# trial_triangular_element = simplex(
# trial_triangular_element.vertices[J[1]],
# trial_triangular_element.vertices[J[2]],
# trial_triangular_element.vertices[J[3]])

# test_sign = Combinatorics.levicivita(I)
# trial_sign = Combinatorics.levicivita(J)

# σ = momintegrals_sign(op, test_sign, trial_sign)

# igd = pulled_back_integrand(op, test_local_space, trial_local_space,
# test_triangular_element, trial_triangular_element)
# G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)
# for j ∈ 1:3, i ∈ 1:3
# out[i,j] += G[K[i],L[j]] * σ
# end

# nothing
# end

# function momintegrals!(op::Helmholtz3DOp,
# test_local_space::LagrangeRefSpace{<:Any,0}, trial_local_space::LagrangeRefSpace{<:Any,1},
# test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

# I, J, K, L = SauterSchwabQuadrature.reorder(
# test_triangular_element.vertices,
# trial_triangular_element.vertices, strat)

# test_triangular_element = simplex(
# test_triangular_element.vertices[I[1]],
# test_triangular_element.vertices[I[2]],
# test_triangular_element.vertices[I[3]])

# trial_triangular_element = simplex(
# trial_triangular_element.vertices[J[1]],
# trial_triangular_element.vertices[J[2]],
# trial_triangular_element.vertices[J[3]])

# test_sign = Combinatorics.levicivita(I)
# trial_sign = Combinatorics.levicivita(J)

# σ = momintegrals_sign(op, test_sign, trial_sign)

# igd = pulled_back_integrand(op, test_local_space, trial_local_space,
# test_triangular_element, trial_triangular_element)
# G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)

# for i ∈ 1:3
# out[1,i] += G[L[i]] * σ
# end

# nothing
# end

# function momintegrals!(op::Helmholtz3DOp,
# test_local_space::LagrangeRefSpace{<:Any,1}, trial_local_space::LagrangeRefSpace{<:Any,0},
# test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

# I, J, K, L = SauterSchwabQuadrature.reorder(
# test_triangular_element.vertices,
# trial_triangular_element.vertices, strat)

# test_triangular_element = simplex(
# test_triangular_element.vertices[I[1]],
# test_triangular_element.vertices[I[2]],
# test_triangular_element.vertices[I[3]])

# trial_triangular_element = simplex(
# trial_triangular_element.vertices[J[1]],
# trial_triangular_element.vertices[J[2]],
# trial_triangular_element.vertices[J[3]])

# test_sign = Combinatorics.levicivita(I)
# trial_sign = Combinatorics.levicivita(J)

# σ = momintegrals_sign(op, test_sign, trial_sign)

# igd = pulled_back_integrand(op, test_local_space, trial_local_space,
# test_triangular_element, trial_triangular_element)
# G = SauterSchwabQuadrature.sauterschwab_parameterized(igd, strat)

# for i ∈ 1:3
# out[i,1] += G[K[i]] * σ
# end

# nothing
# end

# function momintegrals_sign(op::HH3DSingleLayerFDBIO, test_sign, trial_sign)
# return 1
# end
# function momintegrals_sign(op::HH3DDoubleLayerFDBIO, test_sign, trial_sign)
# return trial_sign
# end
# function momintegrals_sign(op::HH3DDoubleLayerTransposedFDBIO, test_sign, trial_sign)
# return test_sign
# end
# function momintegrals_sign(op::HH3DHyperSingularFDBIO, test_sign, trial_sign)
# return test_sign * trial_sign
# end
Loading

0 comments on commit 0caed85

Please sign in to comment.