# Couple equations on two meshes via DG penalty terms

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

``````mesh1 = Make1DMesh(ndofs)
fes1 = L2(mesh1, order=order, dgjumps=True)
fes = fes1*fes1

u,v = fes.TnT()

diffusion += -u[1]*v[0] * ds("right", skeleton=True)
``````

However, this code evaluates both u[1] and v[0] at the rhs and thus creates an matrix entry at the wrong position.

My intention would be to have u[1] evaluated at the right boundary of the first mesh and v[0] on the left of the second, i.e. something like

``````diffusion += -u[1]("left")*v[0]("right") * ds(skeleton=True)
``````

which should result in an entry if the “off-diagonal” blocks of the matrix.

I there any way to implement this without manually modifying the resulting matrices?

Many thanks
Jan

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

Hi Christopher,
many thanks. I will have a look and see if this lies within my capabilities

Best
Jan

If you get stuck feel free to come back with quesitons

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

that would work python-only

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?

I there some working example available?

Many thanks
Jan

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

fes = fes1fes2fes1 ?

Many thanks
Jan

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:
if k < 2:
for i in range(kn,(k+1)n):