mortar-dn-4-mpi.edp

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· EDP 代码 · 共 231 行

assert(version>=2.23);
if ( mpisize != 4 ) 
  {
    cout << " sorry number of processeur !=4 " << endl;
    exit(1);
  }

//  Mortar  (4 sub domain) 
//  with matrix -et Precon Conjugade Gradient --
//    Neuman -> Dirichlet . 
//  -------------------------------

func f=1+x+y;
real g=1; 
int withprecon=1; 
macro Grad(u) [ dx(u), dy(u) ] //
int nbsd=4;

macro Psd(U) U[0],U[1],U[2],U[3] //

int labext= nbsd+1;
real meshsize=0.03; 
real meshsizem=meshsize*1.5; 
bool noconforme=0;

mesh Tha,Thm,Thmm;
int [int] regsd(4);	


if(mpirank==0)
  {
    
    mesh Thacopy; //hack to dcl Tha ouside of mortar-msh.hdp
    {  
      include "mortar-msh.hdp"  
	Thacopy=Tha;
      regsd=regi;
    }   
    Tha=Thacopy;
    cout << "mortar : " << endl;
    Thm=Tha;
    Thm=adaptmesh(Thm,meshsizem,IsMetric=1,thetamax=60);
    Thm=emptymesh(Thm);       Thmm=Thm;
    Thmm=trunc(Thm,split=4,1); // for fine integration
    Thmm=emptymesh(Thmm);
    plot(Thm,wait=0,ps="mortar-Thm.eps");
  }

// send meshes 
verbosity=1;
broadcast(processor(0),Tha);
broadcast(processor(0),Thm);
broadcast(processor(0),Thmm);
broadcast(processor(0),regsd);

// build sub domaine mesh :
mesh Thi=trunc(Tha,region==regsd[mpirank],split=1);// les sous domaines

if(noconforme)
 Thi=adaptmesh(Thi,meshsize*(1+0.1*mpirank),IsMetric=1,nbvx=100000,thetamax=60);// les sous domaines


fespace Lh(Thm,P1);
fespace RTh(Thm,[P0edge,P0edge]);
 RTh [Nmx,Nmy]; // ne marche pas car la normal 
//                             n'est definie que une un bord
varf  vNN([ux,uy],[nx,ny]) = int1d(Thm,1)(( nx*N.x + ny*N.y)/lenEdge);
Nmx[]= vNN(0,RTh);

// les joint P0 sur le squelette 
// -----   \int [q] l + \int[p] m 
Lh  lh=0,rhsl=0;

fespace Vhi(Thi,P1);
fespace Ehi(Thi,P0edge);
matrix Asd,Csd,PAsd,PIsd,PJsd;
Vhi usd,vsd,rhssd, pusd,bcsd;
Ehi epssd;
   
real tgv=1e30;



varf cci([l],[u]) = int1d(Thmm,1,qforder=3)(l*u*epssd);
varf vepsi(u,v)= int1d(Thi,1,qforder=10)( (Nmx*N.x + Nmy*N.y)*v/lenEdge);

varf vLapMi([ui],[vi],tgv=tgv) =
			   int2d(Thi)(  Grad(ui)'*Grad(vi)  )
			   // + int1d(Thi,1,qfe=qf1pElump)(alpha*ui*vi)
			   +  int2d(Thi) (f*vi)   +  on(labext,ui=g);  

varf vPLapMi([ui],[vi],tgv=tgv) =
			   int2d(Thi)(  Grad(ui)'*Grad(vi)  )
			   // + int1d(Thi,1,qfe=qf1pElump)(alphap*ui*vi)
			   + on(labext,1,ui=0);
;  

varf  vrhsMi(ui,vi) =   on(labext,ui=g);


usd=0;
vsd=0;
epssd[]= vepsi(0,Ehi);
epssd = -real(epssd <-0.00001) + real(epssd >0.00001);

Csd  = cci(Lh,Vhi);
Asd  = vLapMi(Vhi,Vhi,solver=UMFPACK);
PAsd = vPLapMi(Vhi,Vhi,solver=UMFPACK);
matrix IVL=interpolate(Vhi,Lh,inside=1);
//   v = IVL*l 
varf vonext(u,v)=on(labext,u=1);
varf von1(u,v)=on(1,u=1);
real[int]  onext=vonext(0,Vhi);
real[int]  on1=von1(0,Vhi);
on1= on1 ? 1 : 0;
on1 = onext ? 0 : on1; //  remove df of ext
matrix I1(on1);//   matrix    tgv $i\in Gamma_1 \ Gamma_e $ , 0 otherwise
PIsd=  I1*IVL;//  remove of line not on $Gamma_1 \ Gamma_e $
// so PIsd*l  =  tgv * Interpole l on $Gamma_1 \ Gamma_e $
I1.diag=on1;
matrix AA=I1*Asd;//  remove line not on lab 1 
PJsd= IVL'*AA;
rhssd[]=vLapMi(0,Vhi);

varf  vML(u,v) = int2d(Thm)(u*v*1e-10)+int1d(Thm,1)(u*v);
matrix ML=vML(Lh,Lh);


lh[]=0;
int itera=0;

varf  vbc(u,v) =   int1d(Thm,labext)(v);
real[int] lbc(Lh.ndof),lbc0(Lh.ndof);
lbc=vbc(0,Lh);
lbc = lbc ? 0 : 1 ; 

int what; //  to choose which  funct call 1  SkPb 2:PSkPb, 3:end
// the preconditionner 
func real[int] PSkPb(real[int] &l)
{ 
  if(withprecon)
    {	
      if(mpirank==0)  
	{ 
	  what=2;
	  broadcast(processor(0),what); // get from SkBk routine
	} 
      if(what!=2) return l;   	
      
      int verb=verbosity;   verbosity=0;   itera++;
      real[int] ll= ML^-1*l;
      broadcast(processor(0),ll);     
      ll= lbc .* ll; 
      ll *= tgv;
      pusd[] = PAsd^-1*(vsd[]= PIsd* ll);
      ll  = PJsd*pusd[];  
      if(mpirank==0)
	{
	  for (int i=1;i<4;++i)
	    {
	      processor(i) >> l;
	      ll += l; 
	    }
	  l = ML^-1*ll; 	
	  l= lbc .* l; 
	}
      else 
	processor(0) << ll;        
      verbosity=verb; 
    }
  return l ;
};

func real[int] SkPb(real[int] &l)
{ 
  int verb=verbosity;   verbosity=0;   itera++;
  if(mpirank==0 && what!=3) what=1;
  broadcast(processor(0),what);  
  if(what==2) return PSkPb(l); 
  else if (what !=1) return l;
  broadcast(processor(0),l);   	
  
  vsd[]  = rhssd[];
  vsd[] += Csd* l;
  usd[]  = Asd^-1*vsd[];      
  l  = Csd'*usd[]; 	
  l= lbc .* l;

  if(mpirank==0)
   {
    real[int] ll(l.n);
    for (int i=1;i<4;++i)
      {
        processor(i) >> ll;
        l += ll; 
      }
    }
  else 
    processor(0) << l;
  verbosity=verb; 
  return l ;
};





if(mpirank==0)
  {
    verbosity=100;
    lh[]=0;
     LinearCG(SkPb,lh[],eps=1.e-5,nbiter=100,precon=PSkPb);
    what=3;
    SkPb(lh[]);
  }
 else
   while(what!=3)
     SkPb(lh[]);

plot(usd,bb=[[-1,-1],[1,1]],ps="mortar-"+mpirank+".eps");
cout << "Fin CG  " << mpirank <<  endl;

/*
Brochet:examples++-mpi hecht$ (grep -vh showpage  mortar-?.eps;echo showpage) > mortar.eps
Brochet:examples++-mpi hecht$ gv mortar.eps 
*/