Radially symmetric source

General
1

Hi! I am trying to create radially symmetric sources in R^2 for the purpose of building a (more) well-posed inverse problem. My first idea was to create a 1D fes, simulating functions that depend only on the radius r, then extend these to a 2D fes, which I’d use to solve PDEs and so on.

My first simple idea was

from netgen.geom2d import SplineGeometry
from ngsolve.meshes import Make1DMesh

mesh1 = Make1DMesh(10, lambda x: 2*x-1)
geo = SplineGeometry()
geo.AddCircle((0,0),1,bc=“circle”)
mesh2 = Mesh(geo.GenerateMesh(maxh=0.1))

fes1 = L2(mesh1, order=2)
fes2 = L2(mesh2, order=2)

q = GridFunction(fes1)
q.Set(1-x**2)
f = GridFunction(fes2)
f.Set(q)

Well, this obviously fails in that f becomes constant in y, but it has the correct behaviour in x. Is there a way to extend q to a radially symmetric f, or a better way to create a “space” of radially symmetric sources? To be clear, just specifying the sources as radially symmetric CFs will not do, since I will be solving inverse problems into the space of radially symmetric sources.

Best,
Christian

2

Not documented, but if your 1d function is smooth you can use GlobalInterfaceSpace:
something like this:

from netgen.occ import *
from ngsolve import *

r=1
geo = OCCGeometry(Circle((0,0),r).Face(), dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
mesh.Curve(3)

fes1 = comp.GlobalInterfaceSpace(mesh, mapping=sqrt((x**2 + y**2)/r**2), order=10)
fes2 = L2(mesh, order=2)

q = GridFunction(fes1)
r = sqrt(x**2 + y**2)
q.Set(1-r**2)
f = GridFunction(fes2)
f.Set(q)

mapping needs to map your mesh coordinates into parameter space [0,1]

3

Thank you very much for the answer! Very interesting to see such a neat feature.

I am not sure I understand the results I get, though. With your code as-is, it produces a result as below:

image
image981×562 61.7 KB

This does not seem consistent with the function 1-x^2-y^2. However, if I set order = 2 in the 1d fes, I get the below, which does appear to be very close to the exact function 1-x^2-y^2. Is there a clear reason for this?

image
image981×562 65.2 KB

As a side note, if this does work with order = 2, then I have the mapping from L^2 over 1d to radially symmetric sources in L^2 over 2d. The adjoint operation (which I will need) from L^2 over 2d to L^2 over 1d would be something like a polar coordinate integral over angle only. How would one go about to implement that?

4

Sorry my bad, sent you code that I didn’t really test. yes with lower order it looks fine. I’m not yet quite sure why high order is so instable here, will investigate. As mentioned, this space is relatively new and not widely adopted/tested, so handle with care :wink:
globalinterface space creates a global basis on parameterspace [0,1] using the given mapping to parameter space. You can visualize the basis functions of the space

Draw(q, mesh, "q")
for i in range(len(q.vec)):
    q.vec[:] = 0
    q.vec[i] = 1
    Redraw()
    input("next")

you can assemble mass matrix and rhs on the global space and use this to interpolate into the space by solving that system.
works with the symbolic forms:


fes1 = comp.GlobalInterfaceSpace(mesh, mapping=1 - (x**2 + y**2)/r**2, order=3)
u,v = fes1.TnT()
m = BilinearForm(u * v * dx)
m.Assemble()
f = LinearForm(func_on_circle * v * dx)
f.Assemble()
gf = GridFunction(fes1)
gf.vec.data = m.mat.Inverse() * f.vec

gf would then be the average of func_on_circle over angle.