dfft.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

// Example of dynamic function load
// --------------------------------
// $Id: dfft.edp,v 1.4 2005/11/01 16:08:57 hecht Exp $
//   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
 // $ -\Delta 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 ?: \label{?:} 
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);