Issue in the Additive Schwarz method

General Discussion
1

Hi

According to the following Additive Schwarz method code, max (uold1,uold2) is to be set on the internal boundaries of both subproblems. But as you know, we can not compare and take the maximum of two vectors unless all similar components of them are well compared in an appropriate way. For example, we can say that uold1 < uold2 if uold1(i) < uold2(i) for all i components on both boundaries. Since the length of uold1 is different from uold2 because of the different number of nodes on the internal boundaries, the code does not run properly.
My question is, Does FreeFem++ automatically compare between all components of the two vectors uold1 and uold2? If yes, then I am in the safe side. If no, then there is a mistake in the code and the results are not correct?
Here is my code and appreciate your help.

2

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.25uold1^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;
};

3

Dear developers

Any help to my problem please?
I want to identify and print the vector values of ‘uold1’ and ‘uold2’ on the internal boundaries labeled 2 on the first problem Schw1 and 4 on the second problem Schw2, then find the maximum of their component values.

For more understanding, we need to find uold1 values, and interpolate uold2 on the first problem internal boundary with label 2. We then compare between them. On the second problem, we need to find uold2 values, and interpolate uold1 on the internal boundary with label 4, and then compare between them.
Attached is my failed try
Additive Schwarz method.edp (1.9 KB)

Thanks in advance.

4

Sorry, I forgot to explain you the technique of interpolation I use it here. In order to get the nodal values on the first internal boundary labeled 2, I only need to interpolate the nodal values of uold2 of the elements surrounding those first nodal values lying on the boundary. That is, the interpolation operator uses only nodal values of uold2 from the the second triangulation to get the nodal values on the first internal boundary.
Hope this is clear enough and get a help from you accordingly.
Thanks

5

Hi developers
Is there any help to my above detailed issue, please? I am still searching for a solution.
Thank you for all your helpful efforts in advance.