Merge pull request #125 from sbadrian/feature/enable_mixed_Helmholtz_…

…Discretizations

Feature/enable mixed helmholtz discretizations

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

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

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
+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