This is not quite right certainly not optimal but thought it was an intersting

approach. I guess you could define two meshes but then you need to

fabricate the boundary equation probably with some messy interpolation matrix

although that is easy to generate in FF with “interpolate” command. It is probably

possible to just use one dof and a discontinuous element assigning the boundary

to one or the other variables. However, I’ve never used those and myabe

someone can comment. This solution appears to work although it wastes

a lot of dofs. And note the picture just sums the variables together causing the

boundary “value” to be inflated ( it is not a sum ) but ok for this picture.

I made up BC and relation between vars but this appears to work and i simple

to write. Thoughts? Maybe someone can modify it to use just one dof if you

have the same element types, this example used diffeing ones

innerline.edp (1.4 KB)

```
load "medit"
macro buildouter()
border left(t=0, 1){x=-szx/2; y=szy*(t-.5);label=1; }
border top(t=0, 1){x=szx*(t-.5);y=szy/2; label=2; }
border right(t=0, 1){x=szx/2; y=-szy*(t-.5);label=3; }
border bottom(t=0, 1){x=-szx*(t-.5);y=-szy/2; label=4; }
// EOM
real szx=2;
real szy=1;
int nx=20;
int ny=20;
buildouter
border middle(t=0, 1){x=0; y=-szy*(t-.5);label=5; }
mesh Th = buildmesh(left(-ny)+top(-nx)+right(-ny)+bottom(-nx)
+ middle(ny));
//medit("Th",Th);
//plot(Th,wait=1,fill=1);
//medit("Th",Th);
macro Op(xx,vxx) (dx(xx)*dx(vxx)+dy(xx)*dy(vxx)) //
fespace Vh1(Th,P1);
fespace Vh2(Th,P2);
Vh1 a1,va1,b1=0;
Vh2 a2,va2,b2=0;
// https://doc.freefem.org/examples/mesh-generation.html
int rl=Th(-szx/2,0).region;
int rr=Th(szx/2,0).region;;
cout<<" rl="<<rl<<" rr="<<rr<<endl; cout.flush;
problem both([a1,a2],[va1,va2],solver=sparsesolver,tgv=-1)
//problem both([a1,a2],[va1,va2],solver=CG)
= int2d(Th,rl)(Op(a1,va1))+int2d(Th,rl)(b1*va1)
+ int2d(Th,rr)(Op(a2,va2))+int2d(Th,rr)(b2*va2)
// otherwise use CG
+int2d(Th,rr)(a1*va1)+int2d(Th,rl)(a2*va2)
// a1 to a2 boundary
+int1d(Th,5)(a2*va2)-int1d(Th,5)(a1*2.0*va2)
+int1d(Th,5)(a2*va1)-int1d(Th,5)(a1*2.0*va1)
+on(1,a1=2*y)+on(3,a2=-3*y)
;
both;
Vh1 net=a1+a2;
plot(a1,cmm="a1",wait=1,fill=true,value=1);
plot(a2,cmm="a2",wait=1,fill=true,value=1);
plot(net,cmm="net",wait=1,fill=true,value=1);
medit("net",Th,net);
```