schwarz-nm.edp

bool withmetis=0; 
bool RAS=0; 
int sizeoverlaps=2; // size off overlap 
int nn=5,mm=5; 

func bool AddLayers(mesh & Th,real[int] &ssd,int n,real[int] &unssd) 
{ 
  //  build a continuous function  uussd (P1) : 
  //  ssd in the caracteristics function on the input sub domain. 
  //  such that : 
  //   unssd = 1 when   ssd =1; 
  //   add n layer of element (size of the overlap) 
  //   and unssd = 0 ouside of this layer ... 
  // --------------------------------- 
  fespace Vh(Th,P1); 
  fespace Ph(Th,P0); 
  Ph s; 
  assert(ssd.n==Ph.ndof); 
  assert(unssd.n==Vh.ndof); 
  unssd=0; 
  s[]= ssd; 
  //  plot(s,wait=1,fill=1); 
  Vh u; 
  varf vM(u,v)=int2d(Th,qforder=1)(u*v/area); 
  matrix M=vM(Ph,Vh); 

  for(int i=0;i<n;++i) 
    { 
      u[]= M*s[]; 
      // plot(u,wait=1); 
      u = u>.1; 
      // plot(u,wait=1); 
      unssd+= u[]; 
      s[]= M'*u[];//'; 
      s = s >0.1; 
    } 
  unssd /= (n); 
  u[]=unssd; 
  ssd=s[]; 
  return true; 
} 

int withplot=3; 
mesh Th=square(100,100); 
int[int] chlab=[1,1  ,2,1  ,3,1  ,4,1  ]; 
Th=change(Th,refe=chlab); 
int npart= nn*mm; 
fespace Ph(Th,P0); 
fespace Vh(Th,P1); 

Ph  part; 
Vh  sun=0,unssd=0; 
Ph xx=x,yy=y,nupp; 
part = int(xx*nn)*mm + int(yy*mm); 
//plot(part,wait=1); 
if(withmetis) 
  { 
    load "metis"; 
    int[int] nupart(Th.nt); 
    metisdual(nupart,Th,npart); 
    for(int i=0;i<nupart.n;++i) 
      part[][i]=nupart[i]; 
  } 
if(withplot>1) 
plot(part,fill=1,cmm="dual",wait=1); 
mesh[int] aTh(npart); 
mesh Thi=Th; 
fespace Vhi(Thi,P1); 
Vhi[int] au(npart),pun(npart); 
matrix[int] Rih(npart); 
matrix[int] Dih(npart); 
matrix[int] aA(npart); 
Vhi[int] auntgv(npart); 
Vhi[int] rhsi(npart); 

for(int i=0;i<npart;++i) 
  { 
    Ph suppi= abs(part-i)<0.1; 
    AddLayers(Th,suppi[],sizeoverlaps,unssd[]); 
    Thi=aTh[i]=trunc(Th,suppi>0,label=10,split=1); 
    Rih[i]=interpolate(Vhi,Vh,inside=1); //  Vh -> Vhi 
    if(RAS) 
      { 
        suppi= abs(part-i)<0.1; 
        varf vSuppi(u,v)=int2d(Th,qforder=1)(suppi*v/area); 
        unssd[]= vSuppi(0,Vh); 
        unssd = unssd>0.; 
        if(withplot>19) 
          plot(unssd,wait=1); 
      } 
    pun[i][]=Rih[i]*unssd[]; 
    sun[] += Rih[i]'*pun[i][];//'; 
    if(withplot>9) 
      plot(part,aTh[i],fill=1,wait=1); 
  } 
plot(sun,wait=1,dim=3,fill=1); 
for(int i=0;i<npart;++i) 
  { 
    Thi=aTh[i]; 
    pun[i]= pun[i]/sun; 
    if(withplot>8) 
      plot(pun[i],wait=1); 
  } 

//  verif partition of unite 

macro Grad(u) [dx(u),dy(u)]//EOM 
  sun=0; 

for(int i=0;i<npart;++i) 
  { 
    cout << " build part :" << i << "/" << npart << endl; 
    Thi=aTh[i]; 
    varf va(u,v) = 
      int2d(Thi)(Grad(u)'*Grad(v))//') 
      +on(1,u=1) + int2d(Thi)(v) 
      +on(10,u=0) ; 


    aA[i]=va(Vhi,Vhi); 
    set(aA[i],solver=UMFPACK); 
    rhsi[i][]= va(0,Vhi); 
    Dih[i]=pun[i][]; 
    real[int]  un(Vhi.ndof); 
    un=1.; 
    real[int] ui=Dih[i]*un; 
    sun[] += Rih[i]'*ui;;//'; 
    varf vaun(u,v) = on(10,u=1); 
    auntgv[i][]=vaun(0,Vhi); // store arry of tgv on Gamma intern. 
  } 
if(withplot>5) 
  plot(sun,fill=1,wait=1); 
cout << sun[].max << " " << sun[].min<< endl; 
// verification of the partition of the unite. 
assert( 1.-1e-9 <= sun[].min  && 1.+1e-9 >= sun[].max); 

int nitermax=1000; 
{ 
  Vh un=0; 
  for(int iter=0;iter<nitermax;++iter) 
    { 
      real err=0; 
      Vh un1=0; 
      for(int i=0;i<npart;++i) 
        { 
          Thi=aTh[i]; 
          real[int] ui=Rih[i]*un[];//'; 
          //{   Vhi uuu; uuu[]=ui;      plot(uuu,wait=1);} 
          real[int] bi = ui .* auntgv[i][]; 
          bi = auntgv[i][] ? bi :  rhsi[i][]; 
          ui=au[i][]; 
          ui= aA[i] ^-1 * bi; 
          //{   Vhi uuu; uuu[]=ui;      plot(uuu,wait=1);} 
          bi = ui-au[i][]; 
          err += bi'*bi;//'; 
          au[i][]= ui; 
          bi = Dih[i]*ui; 
          un1[] += Rih[i]'*bi;//'; 
        } 
      err= sqrt(err); 
      cout << iter << " Err = " << err << endl; 
      if(err<1e-3) break; 
      //    plot(un1,wait=1); 
      un[]=un1[]; 
      if(withplot>2) 
        plot(au,dim=3,wait=0,cmm=" iter  "+iter,fill=1 ); 
    } 
  plot(un,wait=1,dim=3,fill=1); 
}