ilut_edp.edp

load "ilut" 
mesh Th=square(10,10); 
fespace Qh(Th,P1); 
fespace Vh(Th,P2); 
real nu=1; 
varf mp(p,q)=int2d(Th)(p*q/nu); 

matrix Mp=mp(Qh,Qh); 
real[int] diagMp(Qh.ndof); diagMp=Mp.diag; diagMp=1/nu*diagMp; 

// We are putting this here to workaround a bug in freefem++ 
real[int] ilutOUT(Vh.ndof+Qh.ndof); 

// Uses ILUT for the velocity part and the mass matrix for the pressure part 
func real[int] ilutMp(real[int] & xx) { 
	real[int] uuin(Vh.ndof); 
	real[int] ppin(Qh.ndof); 

	real[int] uuout(Vh.ndof); 
	real[int] ppout(Qh.ndof); 

	[uuin,ppin]=xx; 

	for(int k=0;k<Qh.ndof;++k)  { 
		ppout[k]=ppin[k]/diagMp[k]; 
	} 
	//ppout=Mp^-1*ppin; 
	uuout=applyIlutPrecond(uuin); 

	ilutOUT=[uuout,ppin]; 
	return ilutOUT; 
} 

// Uses ILUT for the velocity and I for the pressure 
func real[int] ilut(real[int] & xx) { 
	ilutOUT=applyIlutPrecond(xx); 
	return ilutOUT; 
} 

// Dummy do-nothing preconditioner 
func real[int] dummy(real[int] & xx) { 
	ilutOUT=xx; 
	return ilutOUT; 
}