// The Schwarz alternating method for a nonlinear elliptic problem on two subdomains //where the original domain is [0,1]*[0,1]. The forcing term is f(u)= -u/(1+0.25*u^2) int i=0; //zero iteration func g= x+y; //DIRICHELET boundary condition func c=1; //The function c(x), coeficent of u // Number of vertices array for the 4 meshes of both subdomains real[int] NDL1(4), NDL2(4); real norm1, norm2; real Globalnorm; //real uh; ///////////////////////////////////////////////////////////////////////////// //SUB-DOMAIN1 borders border AB1(t=0, 5./8.){x=t; y=0; label=1;} border BC1(t=0, 1){x=5./8.; y=t; label=2;} border CD1(t=5./8., 0){x=t; y=1; label=3;} border DA1(t=1, 0){x=0; y=t; label=4;} //SUB-DOMAIN2 borders border AB2(t=3./8., 1){x=t; y=0; label=1;} border BC2(t=0, 1){x=1; y=t; label=2;} border CD2(t=1, 3./8.){x=t; y=1; label=3;} border DA2(t=1, 0){x=3./8.; y=t; label=4;} //4 MESHES OF BOTH SUB-DOMAINS for (int n=0; n<=3; n++){ mesh Th1=buildmesh(AB1(10*2^n)+BC1(15*2^n)+CD1(10*2^n)+DA1(15*2^n)); mesh Th2=buildmesh(AB2(7*2^n)+BC2(10*2^n)+CD2(7*2^n)+DA2(10*2^n)); //Finite element spaces fespace Vh1(Th1, P1); Vh1 uh1, vh1, f1, uold1, diff1; //To get the size of the maximum triangle of Th1 fespace Ph1(Th1,P0); Ph1 h1=hTriangle; fespace Vh2(Th2, P1); Vh2 uh2, vh2, f2, uold2, diff2; //To get the size of the maximum triangle of Th1 fespace Ph2(Th2,P0); Ph2 h2=hTriangle; // plot to see the 2 meshes. plot(Th1, Th2, wait=true); //The two Problems f1=-uold1/(1+0.25*uold1^2); problem Schw1(uh1, vh1, init=i, solver=CG, eps=1.e-6) =int2d(Th1)(dx(uh1)*dx(vh1)+dy(uh1)*dy(vh1)) //bilinear term +int2d(Th1)(c*uh1*vh1) +int2d(Th1)(f1*vh1) //right hand side(nonlinear term ) +on(1,3,4, uh1=g)+on(2, uh1=max(uold2,uold1)); //external+internal boundaries f2=-uold2/(1+0.25*uold2^2); problem Schw2(uh2, vh2, init=i, solver=CG, eps=1.e-6) =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(uold2,uold1)); //external+internal boundaries /* uold1=uh1; uold2=uh2; if (uh1<=uh2) uh=uh2; else uh=uh1; */ //Calculation loop for each mesh on both subdomains for (i=0; i<=15; i++){ //solve both problems uold1=uh1; uold2=uh2; Schw1; Schw2; diff1=abs(uh1-uold1); norm1=diff1[].max; diff2=abs(uh2-uold2); norm2=diff2[].max; Globalnorm=max(norm1, norm2); if (Globalnorm<.00001) break; plot(uh1,uh2, wait=true); }; };