# Metric graphs

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):

# second and third edge

for i in range(N + 1):

# add three edges to the mesh

for i in range(N):
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?

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?

using spaces on boundaries/edges should be possible

best Christopher

Hi Christopher,

``````#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.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

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

best, Christopher

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

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