From 95789fa4a7bdd232d7a50f4ecb53ffb5148a6cb9 Mon Sep 17 00:00:00 2001 From: bmergl Date: Thu, 24 Oct 2024 18:39:40 +0200 Subject: [PATCH] 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)