dfft_edp.edp

// Example of dynamic function load 
// -------------------------------- 
// Id 
//   Discret Fast Fourier Transform 
// ------------------------------- 
 load "dfft" 

int nx=32,ny=16,N=nx*ny; 
// warning the fourier space is not exactly the unite square due to periodic condition 
mesh Th=square(nx-1,ny-1,[(nx-1)*x/nx,(ny-1)*y/ny]); 
// warring  the numbering is of the vertices (x,y) is 
// given by   i = x/nx + nx* y/ny 

fespace Vh(Th,P1); 

func f1 = cos(2*x*2*pi)*cos(3*y*2*pi); 
Vh<complex> u=f1,v; 
Vh w=f1; 


Vh  ur,ui; 
//  in dfft the matrix n,m is in row-major order ann array n,m is 
// store j + m* i ( the transpose of the square numbering ) 
 v[]=dfft(u[],ny,-1); 
 u[]=dfft(v[],ny,+1); 
 cout << " ||u||_\infty " << u[].linfty << endl; 
 u[] *= 1./N; // remark: operator /= bug  before version 2.0-3 (FH) 
 cout << " ||u||_\infty " << u[].linfty <<  endl; 
  ur=real(u); 
  plot(w,wait=1,value=1,cmm="w"); 
  plot(ur,wait=1,value=1,cmm="u"); 
 v = w-u; 
cout << " diff = "<< v[].max << " " <<  v[].min << endl; 
assert( norm(v[].max) < 1e-10 &&  norm(v[].min) < 1e-10) ; 
 // -------  a more hard example ----\hfilll 
 // Lapacien en FFT \hfilll 
 //  -Δu = f with biperiodic condition \hfilll 
func f = cos(3*2*pi*x)*cos(2*2*pi*y); // 
func ue =  +(1./(square(2*pi)*13.))*cos(3*2*pi*x)*cos(2*2*pi*y);  // the exact solution 
Vh<complex> ff = f; 
Vh<complex> fhat; 
fhat[] = dfft(ff[],ny,-1); 

Vh<complex> wij; 
// warning in fact we take mode between -nx/2, nx/2 and -ny/2,ny/2 
//  thank to the operator ?:  
wij = square(2.*pi)*(square(( x<0.5?x*nx:(x-1)*nx)) + square((y<0.5?y*ny:(y-1)*ny))); 
wij[][0] = 1e-5; // to remove div / 0 
fhat[] = fhat[]./ wij[];  // 
u[]=dfft(fhat[],ny,1); 
u[] /= complex(N); 
ur = real(u); // the solution 
w = real(ue); // the exact solution 
plot(w,ur,value=1 ,cmm=" ue   ", wait=1); 
w[] -= ur[]; // array sub 
real err= abs(w[].max)+abs(w[].min) ; 
cout << " err = " << err << endl; 
assert( err  < 1e-6);