% sample code of 3D fast Poisson solver (c.f. POICALC)<br>
% written by Weiguo Gao (wggao@fudan.edu.cn)<br>
<br>
% number of grids on each direction<br>
% nz is supposed to be the largest<br>
nx=300; ny=300; nz=300;<br>
% step length<br>
hx=2.*pi/(nx+1); hy=2.*pi/(ny+1); hz=2.*pi/(nz+1);<br>
% square of step length<br>
hx2=hx*hx; hy2=hy*hy; hz2=hz*hz;<br>
<br>
% example from Edouard Oudet, University of Savoie<br>
x=zeros(nx,ny,nz);<br>
for k=1:nz<br>
&nbsp;&nbsp;&nbsp;&nbsp;% right hand side: f(x,y,z)=sin(2*x)*cos(z)<br>
&nbsp;&nbsp;&nbsp;&nbsp;x(:,:,k)=sin(2.*(1:nx)'*hx)*cos(k*hz)*ones(1,ny);<br>
end<br>
% boundary condition: u(x,y,z)=sin(2*x)*cos(z)/5<br>
for k=1:nz, for i=1:nx<br>
&nbsp;&nbsp;&nbsp;&nbsp;x(i,1,k)=x(i,1,k)+sin(2.*i*hx)*cos(k*hz)/(5.*hy2);<br>
&nbsp;&nbsp;&nbsp;&nbsp;x(i,ny,k)=x(i,ny,k)+sin(2.*i*hx)*cos(k*hz)/(5.*hy2);<br>
end, end<br>
x(:,:,1)=x(:,:,1)+sin(2.*(1:nx)'*hx)/(5.*hz2)*ones(1,ny);<br>
x(:,:,nz)=x(:,:,nz)+sin(2.*(1:nx)'*hx)/(5.*hz2)*ones(1,ny);<br>
<br>
dx=(4./hx2)*sin(pi/(2*nx+2.)*[1:nx]).^2;<br>
dy=(4./hy2)*sin(pi/(2*ny+2.)*[1:ny]).^2;<br>
<br>
i=[1:nz 2:nz 1:nz-1]'; j=[1:nz 1:nz-1 2:nz]';<br>
trid=sparse(i,j,-ones(3*nz-2,1)/hz2);<br>
<br>
fprintf('3D fast Poisson solver (Fourier) ... \t');<br>
v=x;<br>
t1=cputime();<br>
for k=1:nz, v(:,:,k)=idst(v(:,:,k)); end<br>
for k=1:nz, v(:,:,k)=idst(v(:,:,k)')'; end<br>
for j=1:ny, for i=1:nx<br>
&nbsp;&nbsp;&nbsp;&nbsp;trid(1:nz+1:nz*nz)=dx(i)+dy(j)+2./hz2;<br>
&nbsp;&nbsp;&nbsp;&nbsp;v(i,j,:)=trid\reshape(v(i,j,:),nz,1);<br>
end, end<br>
for k=1:nz, v(:,:,k)=dst(v(:,:,k)')'; end<br>
for k=1:nz, v(:,:,k)=dst(v(:,:,k)); end<br>
fprintf('%f(s)\n',cputime()-t1);<br>
clear i j k t1 dx dy trid<br>