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

Best
Jan

If you get stuck feel free to come back with quesitons :slight_smile:

you can start with this: GitHub - NGSolve/ngsolve-addon-template: An Minimal NGSolve Addon with Python Bindings

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