From 43903224367ffd7a5119e9b4997456e642888e61 Mon Sep 17 00:00:00 2001 From: Simon Adrian Date: Sun, 24 Mar 2024 22:19:27 +0100 Subject: [PATCH 1/3] Correct momintegration on mixed dual-primal test and trial combinations Currently, mixed discretizations, at least with the scalar primal and dual functions for Helmholtz-type operators are not working due to wrong hardcoded loop boundaries. This is now fixed by using the number of shape functions. --- src/quadrature/sauterschwabints.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/quadrature/sauterschwabints.jl b/src/quadrature/sauterschwabints.jl index b14ad0d6..88e6972f 100644 --- a/src/quadrature/sauterschwabints.jl +++ b/src/quadrature/sauterschwabints.jl @@ -144,12 +144,13 @@ function momintegrals_test_refines_trial!(out, op, Q = restrict(trial_local_space, trial_chart, chart) zlocal = zero(out) + momintegrals!(op, test_local_space, trial_local_space, test_chart, chart, zlocal, qr) - for j in 1:3 - for i in 1:3 - for k in 1:3 + for j in 1:numfunctions(trial_local_space) + for i in 1:numfunctions(test_local_space) + for k in 1:size(Q, 2) out[i,j] += zlocal[i,k] * Q[j,k] end end end end end @@ -194,9 +195,9 @@ function momintegrals_trial_refines_test!(out, op, momintegrals!(op, test_local_space, trial_local_space, chart, trial_chart, zlocal, qr) - for j in 1:3 - for i in 1:3 - for k in 1:3 + for j in 1:numfunctions(trial_local_space) + for i in 1:numfunctions(test_local_space) + for k in 1:size(Q, 2) out[i,j] += Q[i,k] * zlocal[k,j] end end end end From aa4460d065b61d612838e631e56bf5316faf23e1 Mon Sep 17 00:00:00 2001 From: Simon Adrian Date: Mon, 25 Mar 2024 11:20:53 +0100 Subject: [PATCH 2/3] Allow interpolatory dual piecewise constant basis functions These can be useful for mixed discretization of Helmholtz3D type of integral equations. --- src/bases/lagrange.jl | 40 +++++++++++++++++++++++++++++----------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/src/bases/lagrange.jl b/src/bases/lagrange.jl index 5d27f416..271caa1a 100644 --- a/src/bases/lagrange.jl +++ b/src/bases/lagrange.jl @@ -231,15 +231,24 @@ end """ duallagrangecxd0(mesh, jct) -> basis -Build dual Lagrange piecewise constant elements. Boundary nodes are only considered if they are in the interior of `jct`. +Build dual Lagrange piecewise constant elements. Boundary nodes are only considered if +they are in the interior of `jct`. + +The default dual function (`interpolatory=false`) is similar to the one depicted +in Figure 3 of Buffa et al (doi: 10.1090/S0025-5718-07-01965-5), with the +difference that each individual shape function is normalized with respect to +the area so that overall the integral over the dual function is one. + +When `interpolatory=true` is used, the function value is one on the support, and thus, +it gives rise to a partition of unity. """ -function duallagrangecxd0(mesh, jct=CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1)) +function duallagrangecxd0(mesh, jct=CompScienceMeshes.mesh(coordtype(mesh), dimension(mesh)-1); interpolatory=false) vertexlist = interior_and_junction_vertices(mesh, jct) - duallagrangecxd0(mesh, vertexlist) + duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory) end -function duallagrangecxd0(mesh, vertexlist::Vector{Int}) +function duallagrangecxd0(mesh, vertexlist::Vector{Int}; interpolatory=false) T = coordtype(mesh) @@ -252,7 +261,7 @@ function duallagrangecxd0(mesh, vertexlist::Vector{Int}) for (k,v) in enumerate(vertexlist) n = vton[v] F = vtoc[v,1:n] - fns[k] = singleduallagd0(fine, F, v) + fns[k] = singleduallagd0(fine, F, v, interpolatory=interpolatory) push!(pos, verts[v]) end @@ -261,26 +270,35 @@ function duallagrangecxd0(mesh, vertexlist::Vector{Int}) end -function duallagrangecxd0(mesh, vertices::CompScienceMeshes.AbstractMesh{U,1}) where {U} +function duallagrangecxd0(mesh, vertices::CompScienceMeshes.AbstractMesh{U,1}; interpolatory=false) where {U} # vertexlist = Int[v[1] for v in vertices] vertexlist =Int[CompScienceMeshes.indices(vertices, v)[1] for v in vertices] - return duallagrangecxd0(mesh, vertexlist) + return duallagrangecxd0(mesh, vertexlist; interpolatory=interpolatory) end """ - singleduallagd0(fine, F, v) + singleduallagd0(fine, F, v; interpolatory=false) + +Build a single dual constant Lagrange element a mesh `fine`. `F` contains the indices +to cells in the support and v is the index in the vertex list of the defining vertex. + +The default dual function (`interpolatory=false`) is similar to the one depicted +in Figure 3 of Buffa et al (doi: 10.1090/S0025-5718-07-01965-5), with the +difference that each individual shape function is normalized with respect to +the area so that overall the integral over the dual function is one. -Build a single dual constant Lagrange element a mesh `fine`. `F` contains the indices to cells in the support and v is the index in the vertex list of the defining vertex. +When `interpolatory=true` is used, the function value is one on the support, and thus, +it gives rise to a partition of unity. """ -function singleduallagd0(fine, F, v) +function singleduallagd0(fine, F, v; interpolatory=false) T = coordtype(fine) fn = Shape{T}[] for cellid in F # cell = cells(fine)[cellid] ptch = chart(fine, cellid) - coeff = 1 / volume(ptch) / length(F) + coeff = interpolatory ? T(1.0) : 1 / volume(ptch) / length(F) refid = 1 push!(fn, Shape(cellid, refid, coeff)) end From b914087b7a7f8946849b17be0f546da11b7db456 Mon Sep 17 00:00:00 2001 From: Simon Adrian Date: Mon, 25 Mar 2024 11:24:03 +0100 Subject: [PATCH 3/3] Add unit test for mixed discretizations Overall, the observed error is comparable to the standard discretizations. Note, however, that it might be useful (TODO) adding unit tests for individual matrix element assemblies. --- test/test_hh3d_nearfield.jl | 149 +++++++++++++++++++++++++++++++++++- 1 file changed, 148 insertions(+), 1 deletion(-) diff --git a/test/test_hh3d_nearfield.jl b/test/test_hh3d_nearfield.jl index 7cdb04ea..6081d7c7 100644 --- a/test/test_hh3d_nearfield.jl +++ b/test/test_hh3d_nearfield.jl @@ -155,4 +155,151 @@ using Test @test err_EDPDL_field < 0.03 @test err_ENPSL_field < 0.01 @test err_ENPDL_field < 0.02 -end \ No newline at end of file +end + +## Test here some of the mixed discretizatinos +#@testset "Helmholtz potential operators: mixed discretizations with duallagrangecxd0" begin + r = 10.0 + λ = 20 * r + k = 2 * π / λ + + sphere = meshsphere(r, 0.2 * r) + X0 = lagrangecxd0(sphere) + X1 = lagrangec0d1(sphere) + Y0 = duallagrangecxd0(sphere; interpolatory=true) + + S = Helmholtz3D.singlelayer(; gamma=im * k) + D = Helmholtz3D.doublelayer(; gamma=im * k) + Dt = Helmholtz3D.doublelayer_transposed(; gamma=im * k) + N = Helmholtz3D.hypersingular(; gamma=im * k) + + q = 100.0 + ϵ = 1.0 + + # Interior problem + # Formulations from Sauter and Schwab, Boundary Element Methods(2011), Chapter 3.4.1.1 + pos1 = SVector(r * 1.5, 0.0, 0.0) # positioning of point charges + pos2 = SVector(-r * 1.5, 0.0, 0.0) + + charge1 = Helmholtz3D.monopole(position=pos1, amplitude=q/(4*π*ϵ), wavenumber=k) + charge2 = Helmholtz3D.monopole(position=pos2, amplitude=-q/(4*π*ϵ), wavenumber=k) + + # Potential of point charges + + Φ_inc(x) = charge1(x) + charge2(x) + + gD1 = assemble(DirichletTrace(charge1), Y0) + assemble(DirichletTrace(charge2), Y0) + gN = assemble(∂n(charge1), X1) + assemble(BEAST.n ⋅ grad(charge2), X1) + + G = assemble(Identity(), X1, Y0) + o = ones(numfunctions(X1)) + + # Interior Dirichlet problem - compare Sauter & Schwab eqs. 3.81 + M_IDPDL = (-1 / 2 * assemble(Identity(), Y0, X1) + assemble(D, Y0, X1)) # Double layer (DL) + + # Interior Neumann problem + # Neumann derivative from SL potential with deflected nullspace + M_INPSL = (1 / 2 * assemble(Identity(), X1, Y0) + assemble(Dt, X1, Y0)) + G * o * o' * G + + ρ_IDPDL = M_IDPDL \ (-gD1) + ρ_INPSL = M_INPSL \ (-gN) + + pts = meshsphere(0.8 * r, 0.8 * 0.6 * r).vertices # sphere inside on which the potential and field are evaluated + + pot_IDPDL = potential(HH3DDoubleLayerNear(im * k), pts, ρ_IDPDL, X1; type=ComplexF64) + + pot_INPSL = potential(HH3DSingleLayerNear(im * k), pts, ρ_INPSL, Y0; type=ComplexF64) + + # Total field inside should be zero + err_IDPDL_pot = norm(pot_IDPDL + Φ_inc.(pts)) / norm(Φ_inc.(pts)) + err_INPSL_pot = norm(pot_INPSL + Φ_inc.(pts)) / norm(Φ_inc.(pts)) + + # Efield(x) = - grad Φ_inc(x) + Efield(x) = -grad(charge1)(x) + -grad(charge2)(x) + + field_IDPDL = -potential(HH3DHyperSingularNear(im * k), pts, ρ_IDPDL, X1) + field_INPSL = -potential(HH3DDoubleLayerTransposedNear(im * k), pts, ρ_INPSL, Y0) + + err_IDPDL_field = norm(field_IDPDL + Efield.(pts)) / norm(Efield.(pts)) + err_INPSL_field = norm(field_INPSL + Efield.(pts)) / norm(Efield.(pts)) + + # errors of interior problems + @test err_IDPDL_pot < 0.005 + @test err_INPSL_pot < 0.002 + + @test err_IDPDL_field < 0.008 + @test err_INPSL_field < 0.002 +#end + +## Test here some of the mixed discretizatinos +#@testset "Helmholtz potential operators: mixed discretizations with duallagrangec0d1" begin + r = 10.0 + λ = 20 * r + k = 2 * π / λ + + sphere = meshsphere(r, 0.2 * r) + X0 = lagrangecxd0(sphere) + Y1 = duallagrangec0d1(sphere) + + S = Helmholtz3D.singlelayer(; gamma=im * k) + D = Helmholtz3D.doublelayer(; gamma=im * k) + Dt = Helmholtz3D.doublelayer_transposed(; gamma=im * k) + N = Helmholtz3D.hypersingular(; gamma=im * k) + + q = 100.0 + ϵ = 1.0 + + # Interior problem + # Formulations from Sauter and Schwab, Boundary Element Methods(2011), Chapter 3.4.1.1 + pos1 = SVector(r * 1.5, 0.0, 0.0) # positioning of point charges + pos2 = SVector(-r * 1.5, 0.0, 0.0) + + charge1 = Helmholtz3D.monopole(position=pos1, amplitude=q/(4*π*ϵ), wavenumber=k) + charge2 = Helmholtz3D.monopole(position=pos2, amplitude=-q/(4*π*ϵ), wavenumber=k) + + # Potential of point charges + + Φ_inc(x) = charge1(x) + charge2(x) + + gD1 = assemble(DirichletTrace(charge1), X0) + assemble(DirichletTrace(charge2), X0) + gN = assemble(∂n(charge1), Y1) + assemble(BEAST.n ⋅ grad(charge2), Y1) + + G = assemble(Identity(), Y1, X0) + o = ones(numfunctions(X0)) + + # Interior Dirichlet problem - compare Sauter & Schwab eqs. 3.81 + M_IDPDL = (-1 / 2 * assemble(Identity(), X0, Y1) + assemble(D, X0, Y1)) # Double layer (DL) + + # Interior Neumann problem + # Neumann derivative from SL potential with deflected nullspace + M_INPSL = (1 / 2 * assemble(Identity(), Y1, X0) + assemble(Dt, Y1, X0)) + G * o * o' * G + + ρ_IDPDL = M_IDPDL \ (-gD1) + ρ_INPSL = M_INPSL \ (-gN) + + pts = meshsphere(0.8 * r, 0.8 * 0.6 * r).vertices # sphere inside on which the potential and field are evaluated + + pot_IDPDL = potential(HH3DDoubleLayerNear(im * k), pts, ρ_IDPDL, Y1; type=ComplexF64) + + pot_INPSL = potential(HH3DSingleLayerNear(im * k), pts, ρ_INPSL, X0; type=ComplexF64) + + # Total field inside should be zero + err_IDPDL_pot = norm(pot_IDPDL + Φ_inc.(pts)) / norm(Φ_inc.(pts)) + err_INPSL_pot = norm(pot_INPSL + Φ_inc.(pts)) / norm(Φ_inc.(pts)) + + # Efield(x) = - grad Φ_inc(x) + Efield(x) = -grad(charge1)(x) + -grad(charge2)(x) + + field_IDPDL = -potential(HH3DHyperSingularNear(im * k), pts, ρ_IDPDL, Y1) + field_INPSL = -potential(HH3DDoubleLayerTransposedNear(im * k), pts, ρ_INPSL, X0) + + err_IDPDL_field = norm(field_IDPDL + Efield.(pts)) / norm(Efield.(pts)) + err_INPSL_field = norm(field_INPSL + Efield.(pts)) / norm(Efield.(pts)) + + # errors of interior problems + @test err_IDPDL_pot < 0.02 + @test err_INPSL_pot < 0.025 + + @test err_IDPDL_field < 0.02 + @test err_INPSL_field < 0.025 +#end \ No newline at end of file