C++ Girl
  1. C++ Girl
  3. Équation de Poisson utilisant FFTW avec un domaine rectangular
Intereting Posts
Comportement de liste chaînée C ++ SendInput () ne correspond pas à une pression manuelle sur une touche du clavier en C ++? Réduire sur un tableau dans OpenMP Comment tracer le spectre avec FFTW3 / QWT? Envelopper une fonction Python à l’aide de tableaux NumPy en C ++ Essayer de comprendre «pointeur vers membre» Conversion du vecteur STL de la sous-classe en vecteur de la classe de base Convertir une chaîne hexadécimale commençant par «0x» en signature signée en C ++ Est-il sécuritaire de se débarrasser des matières volatiles? Comment trouver quels éléments se trouvent dans le sac en utilisant l’algorithme Knapsack ? Comment savoir sur quel type d’object un pointeur pointe en C ++? Que signifie `classe HelloWorld: public Gtk :: Window`? Comment obtenir une info-bulle multiligne dans MFC Comment cloner un object en C ++? Ou y a-t-il une autre solution? la sémantique de référence c n’est pas explicite dans le code de l’appelant

Équation de Poisson utilisant FFTW avec un domaine rectangular

J’essaie de résoudre l’équation de Poisson en utilisant une bibliothèque FFTW avec un domaine rectangular (-4 <= x <= 4 et -2 <= y <= 2). J'ai le résultat correct si le domaine est carré et il est faux si le domaine est rectangulaire. S'il vous plaît donnez-moi une suggestion. Merci beaucoup. Voici mon code. 

#include  #include  #include  #include  #include  #include  using namespace std; int main() { int N1=64; int N2=32; double pi = 3.141592653589793; double L1 = 8.0; double dx = L1/(double)(N1-1); double L2= 4.0; double dy=L2/(double)(N2-1); std::vector in1(N1*N2,0.0); std::vector in2(N1*N2,0.0); std::vector out1(N1*N2,0.0); std::vector out2(N1*N2,0.0); std::vector X(N1,0.0); std::vector Y(N2,0.0); fftw_plan p, q; int i,j; p = fftw_plan_r2r_2d(N1,N2, in1.data(), out1.data(), FFTW_REDFT00, FFTW_REDFT00, FFTW_EXHAUSTIVE); q = fftw_plan_r2r_2d(N1,N2, in2.data(), out2.data(), FFTW_REDFT00, FFTW_REDFT00, FFTW_EXHAUSTIVE); int l=-1; for(i = 0;i <N1;i++){ X[i] =-4.0+(double)i*dx ; for(j = 0;j<N2;j++){ l=l+1; Y[j] =-2.0+ (double)j*dy ; in1[l]= cos(pi*X[i]) + cos(pi*Y[j]) ; // row major ordering } } fftw_execute(p); l=-1; for ( i = 0; i < N1; i++){ // f = g / ( kx² + ky² ) for( j = 0; j < N2; j++){ l=l+1; double fact=0; in2[l]=0; if(2*i<N1){ fact=((double)i*i); }else{ fact=((double)(N1-i)*(N1-i)); } if(2*j<N2){ fact+=((double)j*j); }else{ fact+=((double)(N2-j)*(N2-j)); } if(fact!=0){ in2[l] = out1[l]/fact; }else{ in2[l] = 0.0; } } } fftw_execute(q); l=-1; double erl1 = 0.; for ( i = 0; i < N1; i++) { for( j = 0; j < N2; j++){ l=l+1; erl1 += fabs( in1[l]- 8.*out2[l]/((double)(N1-1))/((double)(N2-1))); printf("%3d %10.5f %10.5f\n", l, in1[l], 8.*out2[l]/((double)(N1-1))/((double)(N2-1))); } } cout<<"error=" <<erl1 <<endl ; fftw_destroy_plan(p); fftw_destroy_plan(q); fftw_cleanup(); return 0; }

Dans l’espace de Fourier, les fréquences k_x et k_y doivent dépendre de la taille du domaine. Comme le domaine est rectangular, cela devient une question d’importance.

Essayez: 

 double invL1s=1.0/(L1*L1); double invL2s=1.0/(L2*L2); ... if(2*i

La sortie doit être plus proche de celle attendue (sauf pour un facteur d'échelle).