Domain descomposition techniques for interfacial discontinuties

General Discussion
1

Hi to everyone,

I’m following the following example by Mc Bain: https://www.ljll.math.upmc.fr/hecht/ftp/ff++days/2012/McBain.pdf

which works when the flux is the same, this is k_1 \frac{\partial T}{\partial n} = k_2 \frac{\partial T}{\partial n} .

Can we do the same in FreeFem when there is also a jump in the flux? this is

T1 - T2 = h(T1, T2,k_1\frac{\partial T_1}{\partial n},k_2\frac{\partial T_2}{\partial n}) and k_1\frac{\partial T_1}{\partial n} - k_2\frac{\partial T_2}{\partial n} = g(T1,T2)

In particular, i’m trying this code (minimal), but GC doesn’t converge

int j=0;
Fh1 q;
problem Pb1 (T1, S1, init=j)
      = int2d(Th1)(k1*[dx(T1),dy(T1)]'*[dx(S1),dy(S1)])
        -int2d(Th1)(f1*S1)
       +int1d(Th1, inner)(q*S1)
       +on(left, T1=0);

problem Pb2 (T2, S2, init=j)
     =  int2d(Th2)(k1*[dx(T2),dy(T2)]'*[dx(S2),dy(S2)])
        -int2d(Th2)(f2*S2)
       -int1d(Th2 inner)(q*S2 + h*S2)
      +on(right, T2=1);

varf bq(phi, psi, solver=CG) = int1d(Th1, inner)(phi*psi);
matrix B = bq(Fh1, Fh1, solver=CG);

 func real[int] BoundaryProblem (real[int] &l){
    q[] = l; 
    Pb1;
    Pb2;
    j++; 
    S1 = T1 - T2; 
    real[int] q1 = B * S1[];
    S1 = q + k2 * [dx(T2),dy(T2)]'*[N.x,N.y];
    real[int] q2 = B * S1[];
    q[] = q1 - q2;
    return q[];
}

q=0.;
Fh1 variable=0.; 
Pb1;
Pb2;
LinearCG(BoundaryProblem, variable[], nbiter=1e2);
BoundaryProblem(variable[]);

Could you help me please?
Best regards