Projection of SurfaceL2 function onto VectorH1 space.

[Horep]

Suppose I have a SurfaceL2 function g defined on the entire boundary of my mesh.

I now need to solve a Dirichlet problem $$<∇u,∇v>=0, \quad u=g \text{ on } \partial\Omega$$ where g is the boundary condition for the full boundary. The standard way to do this is to decompose $u = u_0 + g_h$, (g_h being some VectorH1 projection of g) and then solve $$<∇u0,∇v>=-<∇g{h},∇v>, u_0=0\text{ on } \partial\Omega,$$ which is simple enough within NGSolve. I am just not sure how to project g into an H1 function. What is the best way to go about doing this?

[Weggler]

note that ngbem implements boundary value problems.

• To see how to solve an interior Dirichlet bvp, have a look at one of the demos Laplace_DtN_direct or Laplace_DtN_indirect.
• The variational formulation in your poste is a finite element formulation -> check out ngsolve demos for this.

[Horep]

Hi,

The problem is that I must construct the SL matrix, which is (from my understanding) expensive.

For full context, I am trying to implement the following algorithm (from https://doi.org/10.1016/j.camwa.2017.11.028)

Point (i) is simple using NGSolve. Point (ii) is where I am using NGBem (modelled after Laplace_DtN_direct), solving the system $$Mg = (K -0.5 M)u_1$$ which requires the DL potential, and I think I have a working piece of code for this, but it yields a SurfaceL2 object g.

Since I don't want to use the SL potential if possible, I would rather project g into VectorH1 so I can use the usual NGSolve method when solving point (iii).

I would otherwise use one of Laplace_DtN_direct or Laplace_DtN_indirect, but these don't seem to take into account when I am using the DL potential as the boundary condition itself (see (49d)). If I was following Laplace_DtN_direct, would this mean I have $$Vu_2 = ( K^2 +0.5 M) u_1,$$ combining both point (ii) and (iii) into a single step?

[Weggler]

yes, I think so: combine point (ii) and (iii) into a single step!

(49c+d) can be solved as in the demos Laplace_DtN provided you know the Dirichlet trace. The paper proposes an explicit representation of the latter in terms of the FEM solution u1 as (K+1/2M) u1. If you have u1, it might work.

[Horep]

Hi,

So, I have u1, it is easily calculated:

def compute_internal(mag_func: GridFunction, mesh: Mesh) -> GridFunction:
    V = H1(mesh, order=1)
    N = NumberSpace(mesh) # used for enforcing zero mean value, constant FE space
    X = V*N
    (u, lam), (v,mu) = X.TnT()
    a = BilinearForm(X)
    a += InnerProduct(Grad(u), Grad(v)) * dx
    a += (u*mu + v*lam) * dx # enforces zero mean value condition
    f = LinearForm(X)
    f += mag_func*Grad(v) * dx

    gfu = GridFunction(X)
    with TaskManager():
        a.Assemble()
        f.Assemble()
        rows,cols,vals = a.mat.COO()
        A = scipy.sparse.csc_matrix((vals,(rows,cols)))
        F = f.vec.FV().NumPy()[:]
        M2 = scipy.sparse.linalg.spilu(A)  # spilu preconditioner
        M = scipy.sparse.linalg.LinearOperator((F.size, F.size), M2.solve)
        gfu.vec.FV().NumPy()[:] = scipy.sparse.linalg.gmres(A, F, M=M)[0]
    return gfu.components[0]

Now, given u1 I attempt to perform the computation of point (ii), but I am still unsure how to properly do it (I need a "mass matrix" on both sides, we are solving $M'g = (K- 0.5M)u_1$), with

def DL_boundary_projection(pot_int: GridFunction, mesh: Mesh) -> GridFunction:
    fesL2 = SurfaceL2(mesh, order=1, dual_mapping=True)
    u, v = fesL2.TnT()
    fesH1 = H1(mesh, order=1)
    uH1, vH1 = fesH1.TnT()
    g = GridFunction(fesL2)
    with TaskManager(): 
        # surface mass matrix
        M_SURF = BilinearForm( u.Trace() * v.Trace() * ds(bonus_intorder=2)).Assemble()
        M_H1 = BilinearForm( uH1 * v.Trace() * ds(bonus_intorder=2)).Assemble()
        K = DoubleLayerPotentialOperator(fesH1, fesL2, intorder=2, leafsize=40, eta=3., eps=1e-4, method="aca")  # <u_H1, v_L2>
        f = ((K.mat - 0.5 * M_H1.mat)*pot_int.vec).Evaluate()

        rows,cols,vals = M_SURF.mat.COO()
        A = scipy.sparse.csc_matrix((vals,(rows,cols)))
        F = f.FV().NumPy()[:]
        g.vec.FV().NumPy()[:] = scipy.sparse.linalg.gmres(A, F)[0]

    return g

Regardless, at this point we now have the function $g$. If we want to continue like we do with NGSolve, we now need to be able to differentiate $g$. Since SurfaceL2 will not have enough regularity for this, I propose taking an H1 function, with zero on the interior nodes, that matches g.Trace() on the boundary. I am unsure of how to perform this step though, hence this post.

[Weggler]

the mapping (K+1/2M)u1 returns a SurfaceL2 trace and it's not suited for a standard BEM formulation. You should consider (iii) FEM formulation. You should rather look into ngsolve demos and ask the ngsolve community for the projection. Unfortunately I am not yet a ngsolve expert. I have no quick answer.

[Horep]

Okay, I will think about this further, and see if I can come up with a solution.

[Horep]

Unsurprisingly, the solution is to just use the Set method, i.e. gfu.Set(g, BND), provided gfu belongs to a space with Dirichlet conditions on the boundary.