/usr/share/doc/freefem++/examples/examples++-load/dfft.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 | // 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
// $ -\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);
|