Hi,
I am trying to couple several (1d) meshes in order to simulate PDEs on a metric graph.

A minimal example consists of two 1d meshes and a jump term that connects DOFs on the right boundary of the first mesh with those on the left of the second

hi, i don’t think this can be done from python right now, but you can create a small cpp extension and utilize the special elements that are used for contact formulations.

look at contact.hpp/cpp in comp dir in ngsolve. you need basically the same but instead of contact pairs your coupling vertices

and as a first step add special elements with hard coded terms that you want using the contact special elements. From there going to symblic terms you can also use code parts from contact.cpp

Ah maybe a different idea: you could also use the globalinterface with dim 0 and mapping form the 2 vertices to origin to glue together your points with a hybrid technique, similar to HDG

you find the space in ngsolve.comp.GlobalInterfaceSpace

Dear Christopher,
many thanks for your effort. I python-only solution would be great, however, it is not so easy for me to understand what the globalinterfacespace actually does, e.g. how would the mapping look like/which format is appropriate?

Dear Christopher,
sorry to bother again, but would it be possible to provide a minimal example for the use of globalinterfacespace? I think I now understand how it should work, but I am not sure how the mapping would be defined on the product space…

I guess the starting point would be
fes1 = L2(mesh1, order=order, dgjumps=True)
fes2 = globalinterfacespace

Ok what I thought works is currently not yet implemented (GlobalInterfaceSpace on BBND elements). So I think the easier way would be to create a cpp extension and add integrators like the ones in the contact.*pp

Dear Christopher,
after meeting André Massing by chance yesterday, he told me that he had a similar isse (yet in 3d) here he simply used HDG directly, which worked out of the box.

I tried the same, assigning, in a point where several edges meed, the same facet DOF to all adjacent elements and it turned out that this works for my setting as well.

For completeness, this is the (admittedly very handcrafted) mesh-generation routine for a y-shaped graph. After that, I just use HDG as in the NgSolve documentation and the assembly routines worked perfectly fine.

n = number of dofs on each edge
n_edges = 3 in this example
for k in range(n_edges):
for i in range(n+1):
x = i/n
if mapping:
(x,y) = mapping(x,k) # Mapping maps to point on respective edge
if k > 0 and i == 0:
print(“Point skipped”)
else:
pids.append (mesh.Add (MeshPoint(Pnt(x, y, 0))))
if k < 2:
for i in range(kn,(k+1)n):
mesh.Add(Element1D([pids[i],pids[i+1]],index=idx_inner))
else:
mesh.Add(Element1D([pids[n],pids[kn+1]],index=idx_inner))
for i in range(kn+1,(k+1)*n):
mesh.Add(Element1D([pids[i],pids[i+1]],index=idx_inner))

mesh.Add (Element0D( pids[0], index=idx_left)) #Add boundary nodes at unconnected ends of each edge
mesh.Add (Element0D( pids[2n], index=idx_right1))
mesh.Add (Element0D( pids[3n], index=idx_right2))