Metric graphs

General
1

Dear all,
I am trying to solve a PDE on a metric graph, i.e. a graph where to each edge a one-dimensional is associated.

To model this, I tried to create a two-dimensional mesh consisting of three edges made up of one-dimenionsal elements Element1D as follow:

import netgen.gui
import ngsolve
from netgen.meshing import *
from netgen.geom2d import unit_square
ngmesh = Mesh(dim=2)

N=5
pnums =

first edge

for i in range(N + 1):
pnums.append(ngmesh.Add(MeshPoint(Pnt(i / N, 0, 0))))

second and third edge

for i in range(N + 1):
pnums.append(ngmesh.Add(MeshPoint(Pnt(1., i/N, 0))))
pnums.append(ngmesh.Add(MeshPoint(Pnt(1., -1.0*i/ N, 0))))

add three edges to the mesh

for i in range(N):
ngmesh.Add(Element1D([pnums[i], pnums[i + 1]], index=1))
ngmesh.Add(Element1D([pnums[N + i + 1], pnums[N+i + 1 + 1]], index=2))
ngmesh.Add(Element1D([pnums[2N+2 + i], pnums[2N + 2 + i + 1]], index=3))

Does this result in a valid mesh and can I then define bilinear forms on this, simular to the case of surface PDEs?

2

Alternatively, it would also work to solve a PDE on the boundary elements of a 2d-mesh, can this be done with an H^1 space defined on the boundary, only?

3

using spaces on boundaries/edges should be possible

best Christopher

4

Hi Christopher,
thanks for your reply. However, the following code

#import netgen.gui
from ngsolve import *
import ngsolve as ng
from netgen.meshing import *
from netgen.geom2d import unit_square
from netgen.geom2d import SplineGeometry

geo = SplineGeometry()

pnts =[(0,0), (1,0), (0.5,1)]

p1,p2,p3 = [geo.AppendPoint(*pnt) for pnt in pnts]

curves = [[["line",p1,p2],"e1"],
          [["line",p2,p3],"e2"],
          [["line",p3,p1],"e3"]]

[geo.Append(c,bc=bc) for c,bc in curves]

mesh = ng.Mesh(geo.GenerateMesh(maxh=0.1))

fesext = H1(mesh, order=1, definedon="e1")
gfue = GridFunction(fesext)

u, v = fesext.TnT()

a = BilinearForm(fesext, symmetric=True)
a += grad(u).Trace()*grad(v).Trace()*ds
# a += grad(u)*grad(v)*dx # yields same result
a.Assemble()

f = LinearForm(fesext)
f += 1*v*ds
f.Assemble()

gfue.vec.data = a.mat.Inverse(fesext.FreeDofs()) * f.vec

print(gfue.vec)

yields a zero vector as outcome, also f.vec is only zeros. Am I overlooking something or do I need to chose a different space?

thx
Jan

5 
fesext = H1(mesh, order=1, definedon=mesh.Boundaries("e1"))

best, Christopher

6

Dear Christopher,
many thanks, that works!

Follow up: Is there an easy way to visualize the results? Previously, I was using

pnts_vals = [ (x,gfue(x)) for x in pnts if ngsmesh.Contains(x)]
pnts,vals = zip(*pnts_vals)

Can I constrain/restrict that to the boundary part “e1” (in my example)? I tried

for el in ngmesh.Elements1D():
    print(ngmesh.GetBCName(el.index-1))
    for vert in el.vertices:
        print("val at " + str(ngmesh[vert]) + " is " + str(gfue(ngmesh[vert][0],ngmesh[vert][1],ngmesh[vert][2])))

but without success

Many thanks
Jan

7

you need to pass the BND flag to get the mesh points on the edges

vals = cf(mesh(*pnts, VOL_or_BND=BND))

i’d do visualization then in matplotlib or plotly