# Domain descomposition techniques for interfacial discontinuties

**URL:** https://community.freefem.org/t/domain-descomposition-techniques-for-interfacial-discontinuties/2216
**Category:** General Discussion
**Created:** [January 12, 2023, 9:15pm UTC](https://community.freefem.org/t/domain-descomposition-techniques-for-interfacial-discontinuties/2216 "2023-01-12T21:15:37Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![RodZ](https://avatars.discourse-cdn.com/v4/letter/r/ad7895/32.png) [@RodZ](https://community.freefem.org/u/RodZ)
#### Post date: [January 12, 2023, 9:15pm UTC](https://community.freefem.org/t/domain-descomposition-techniques-for-interfacial-discontinuties/2216/1 "2023-01-12T21:15:37Z")

</div>

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](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

```auto
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