Couple equations on two meshes via DG penalty terms

General
1

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

2

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

3

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

Best
Jan

4

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

5

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

6

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

7

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

8

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

9

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

10

Ah super. I thought that your meeting points were spacially separated. Yes if they are together you can do it like this.