int n=0; //iteration

int i=0; //number of meshes

func g= x+y; //DIRICHELET boundary condition

func c=1; //The function c(x), coeficent of u

real x0=0, x1=12./16., y0=0, y1=1; //SUB-DOMAIN1 Dimensions

real xx0=4./16., xx1=1, yy0=0, yy1=1; //SUB-DOMAIN2 Dimensions

real norm1, norm2;

real Globalnorm;

//SUB-DOMAIN1 borders

border AB1(t=x0, x1){x=t; y=y0; label=1;};

border BC1(t=y0, y1){x=x1; y=t; label=2;};

border CD1(t=x1, x0){x=t; y=y1; label=3;};

border DA1(t=y1, y0){x=x0; y=t; label=4;};

//SUB-DOMAIN2 borders

border AB2(t=xx0, xx1){x=t; y=yy0; label=1;};

border BC2(t=yy0, yy1){x=xx1; y=t; label=2;};

border CD2(t=xx1, xx0){x=t; y=yy1; label=3;};

border DA2(t=yy1, yy0){x=xx0; y=t; label=4;};

mesh Th1=square(2, 2, [x0+(x1-x0)*x, y0+(y1-y0)*y]); //initial Rectangular mesh of SUB-DOMAIN1

mesh Th2=square(2, 2, [xx0+(xx1-xx0)*x, yy0+(yy1-yy0)*y]); //initial Rectangular mesh of SUB-DOMAIN2

//4 MESHES OF BOTH SUB-DOMAINS (Two choices of meshings)

for (int i=0; i<=3; i++){

Th1 = adaptmesh(Th1,(1./16.)*(1./(2.^i)), IsMetric=1, thetamax=89);*

Th2 = adaptmesh(Th2,(1./12.)(1./(2.^i)), IsMetric=1, thetamax=89);

//Finite element spaces

fespace Vh1(Th1, P1);

Vh1 uh1, vh1, f1, uold1, diff1, err1, uold;

fespace Vh2(Th2, P1);

Vh2 uh2, vh2, f2, uold2, diff2, err2;

// plot to see the 2 meshes.

plot(Th1, Th2, wait=true);

//The two Problems

uold1=0 ; uold2=0;

f1=-uold1/(1+0.25*uold1^2);*

problem Schw1(uh1, vh1, init=i, solver=CG, eps=1.e-8)

=int2d(Th1)(dx(uh1)*dx(vh1)+dy(uh1)**dy(vh1)) //bilinear term*

+int2d(Th1)(cuh1vh1)

+int2d(Th1)(f1vh1) //right hand side(nonlinear term )

+on(1,3,4, uh1=g)+on(2, uh1=max(uold1,uold2)); //external+internal boundaries

```
f2=-uold2/(1+0.25*uold2^2);
problem Schw2(uh2, vh2, init=i, solver=CG, eps=1.e-8)
=int2d(Th2)(dx(uh2)*dx(vh2)+dy(uh2)*dy(vh2)) //bilinear term
+int2d(Th2)(c*uh2*vh2)
+int2d(Th2)(f2*vh2) //right hand side (nonlinear term)
+on(1,2,3, uh2=g)+on(4, uh2=max(uold1,uold2)); //external+internal boundaries
```

//Calculation loop for each mesh on both subdomains

for (int n=0; n<=10; n++){

//solve both problems

uold1=uh1 ; uold2=uh2;

Schw1;

int NbVertices1=Th1.nv;

for (int j = 0; j <= NbVertices1; j++){

int L1 = Th1(j).label;

real x1= Th1(j).x; real y1= Th1(j).y;

if (L1 == 2) {

cout << " At “<< j << " : " << x1 << " , " << y1 << " , " << uold1[j]<< endl;

};

};

Schw2;

int NbVertices2=Th2.nv;

for (int k = 0; k <= NbVertices2; k++){

int L2 = Th2(k).label;

real x2= Th2(k).x; real y2= Th2(k).y;

if (L2 == 4) {

cout << " At “<< k << " : " << x2 << " , " << y2 << " , " << uold2[k]<< endl;

};

};

int l , m ;

l = min(NbVertices1,NbVertices2);

m = max(NbVertices1,NbVertices2);

for (l = min(NbVertices1,NbVertices2); l <= m; l++){

if (NbVertices1 < NbVertices2) {

uold1[l] = 0;}

else{

uold2[l] = 0;}

};

cout << uold1[m] <<” --and-- “<< uold2[m] << endl;

uold = max(uold1[m],uold2[m]);

cout << “uold = max(uold1,uold2) @ iteration”<<”(”<<n<<“) =”<< uold << endl;

diff1=abs(uh1-uold1);

norm1=diff1.max;

diff2=abs(uh2-uold2);

norm2=diff2.max;

Globalnorm=max(norm1, norm2);

if (Globalnorm<1e-8) break;

cout<<“Difference error on subdomain 1 @ iteration”<<“(”<<n<<“) =”<<norm1<<endl;

cout<<“Difference error on subdomain 2 @ iteration”<<“(”<<n<<“) =”<<norm2<<endl;

cout<<“comparison between the approximate and exact solutions at some certain points”<<endl;

};

cout<<“Numerical Results of DISRETIZATION No.”<<i<<endl;

plot(uh1,uh2, wait=true, fill=true);

cout<<" Number of iterations is:"<<n<<endl;

err1=abs(uh1-uh0);

maxerror1[i]=err1.max;

cout <<“GLOBAL maximum error on subdomain 1 is:” << maxerror1[i] <<endl;

err2=abs(uh2-uh0);

maxerror2[i]=err2.max;

cout <<“GLOBAL maximum error on subdomain 2 is:” << maxerror2[i] <<endl;

};