Dear All,
I’m using BEM add-on library of NGSolve. I’m trying to solve
S[f]={\rm Sf} on unit sphere in two different ways (the Laplace single layer potential is H^{-1/2}-elliptic, so I believe this problem is well-posed given {\rm Sf}):
I define the geometry and mesh
geo = CSGeometry()
sphere = Sphere(Pnt(0,0,0), 1)
geo.Add(sphere.col([0,0,0]))
mesh = geo.GenerateMesh(maxh=0.4, perfstepsend=MeshingStep.MESHSURFACE) # no non-free dof
mesh = Mesh(mesh)
mesh.Curve(2)
and function space and SIO
fes = H1(mesh, order=2, definedon=mesh.Boundaries(".*"))
L_SL = SingleLayerPotentialOperator(fes, intorder=5, leafsize=20, eta=3., eps=1e-4, method="aca", testhmatrix=False)
, and then I define the function f:
f = GridFunction(fes)
f.Set(CF(x+y**2), definedon = mesh.Boundaries('.*'))
Way 1:
I use “.GetPotential()” to compute {\rm Sf}
Sf = GridFunction(fes)
Sf.Set(L_SL.GetPotential(f), definedon = mesh.Boundaries('.*'))
Way 2:
I use {\rm Sf} obtained in “Way 1” to compute f by manually solving the H-matrix linear system of equations
m_id = IdentityMatrix(fes.ndof, complex=False)
f_sol = GMRes(A=L_SL.mat, b=Sf.vec, pre=m_id, tol=1e-8, maxsteps=400, printrates=False)
But it turns out f and f\_sol are far from the same. The second approach is widely used in the tutorial, and I also believe the evaluation function in the first approach is accurate.
Thank you in advance!