function cppspm3d
% cppspm3d: Chebyshev pseudospectral Poisson solver with penalty method in a 3D rectangular domain.
% Citation: Tzyy-Leng Horng and Chun-Hao Teng, 2012, "An error minimized pseudospectral penalty direct Poisson solver," 
% Journal of Computational Physics, 231(6):2498¡V2509 
% uxx+uyy+uzz=f(x,y,z); a test case with exact solution ue(x,y)=sin(4*pi*X).*sin(4*pi*Y).*sin(4*pi*Z);
% BC: amx*u(-1,y,z)-bmx*ux(-1,y,z)=gmx(y,z), at x=-1,
%       apx*u(1,y,z)+bpx*ux(1,y,z)=gpx(y,z), at x=1,
%       amy*u(x,-1,z)-bmy*uy(x,-1,z)=gmy(x,z), at y=-1,
%       apy*u(x,1,z)+bpy*uy(x,1,z)=gpy(x,z), at y=1.
%       amz*u(x,y,-1)-bmz*uz(x,y,-1)=gmz(x,y), at z=-1,
%       apz*u(x,y,1)+bpz*uz(x,y,1)=gpz(x,y), at z=1.
% PS: This code is adapted from poisson_penalty_3d_benchmarking_final.m
  close all;
  amx=input('Input amx, default is 1: ');
  if isempty(amx)
      amx=1;
  end
  bmx=input('Input bmx, default is 0: ');
  if isempty(bmx)
      bmx=0;
  end
  apx=input('Input apx, default is 1: ');
  if isempty(apx)
      apx=1;
  end
  bpx=input('Input bpx, default is 0: ');
  if isempty(bpx)
      bpx=0;
  end
  amy=input('Input amy, default is 1: ');
  if isempty(amy)
      amy=1;
  end
  bmy=input('Input bmy, default is 0: ');
  if isempty(bmy)
      bmy=0;
  end
  apy=input('Input apy, default is 1: ');
  if isempty(apy)
      apy=1;
  end
  bpy=input('Input bpy, default is 0: ');
  if isempty(bpy)
      bpy=0;
  end
  amz=input('Input amz, default is 1: ');
  if isempty(amz)
      amz=1;
  end
  bmz=input('Input bmz, default is 0: ');
  if isempty(bmz)
      bmz=0;
  end
  apz=input('Input apz, default is 1: ');
  if isempty(apz)
      apz=1;
  end
  bpz=input('Input bpz, default is 0: ');
  if isempty(bpz)
      bpz=0;
  end
  Nvec=input('Input N in vector, default is [16 20 24 32 40 48 64]: ');
  if isempty(Nvec)
      Nvec=[16 20 24 32 40 48 64];
  end
  errl2vec=[]; errinfvec=[];
  for N=Nvec
  N
  Nx=N; Ny=N; Nz=N;
  cwx=ones(Nx+1,1); cwx([1 end])=2; 
  cwy=ones(Ny+1,1); cwy([1 end])=2; 
  cwz=ones(Ny+1,1); cwz([1 end])=2; 
  cxywmat=cwy*cwx.';
  cwmat=[];
  for i=1:Nx+1
      tmp=cxywmat(:,i)*cwz.';
      cwmat=cat(3,cwmat,tmp);
  end
  % computing errorminized tau's in x direction:
  am=amx; bm=bmx; ap=apx; bp=bpx;
  gp=(ap/(N^2-1)-bp)/(2*N^2)+(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  gm=(ap/(N^2-1)-bp)/(2*N^2)-(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  hp=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)+(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  hm=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)-(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  taumoptx=(-1)^N*(ap^2+2*(ap+bp)^2)/((2*(ap+bp)*(am+bm)-am*ap)*gm+(ap^2+2*(ap+bp)^2)*hm)
  taupoptx=-(am^2+2*(am+bm)^2)/((am^2+2*(am+bm)^2)*gp+(2*(am+bm)*(ap+bp)-am*ap)*hp)
  % computing errorminized tau's in y direction:
  am=amy; bm=bmy; ap=apy; bp=bpy;
  gp=(ap/(N^2-1)-bp)/(2*N^2)+(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  gm=(ap/(N^2-1)-bp)/(2*N^2)-(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  hp=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)+(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  hm=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)-(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  taumopty=(-1)^N*(ap^2+2*(ap+bp)^2)/((2*(ap+bp)*(am+bm)-am*ap)*gm+(ap^2+2*(ap+bp)^2)*hm)
  taupopty=-(am^2+2*(am+bm)^2)/((am^2+2*(am+bm)^2)*gp+(2*(am+bm)*(ap+bp)-am*ap)*hp)
  % computing errorminized tau's in z direction:
  am=amz; bm=bmz; ap=apz; bp=bpz;
  gp=(ap/(N^2-1)-bp)/(2*N^2)+(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  gm=(ap/(N^2-1)-bp)/(2*N^2)-(ap*(N^2+2)/(N^2*(N^2-4))-bp)/(2*(N^2-1));
  hp=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)+(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  hm=(-1)^(N-1)*(am/(N^2-1)-bm)/(2*N^2)-(-1)^N*(am*(N^2+2)/(N^2*(N^2-4))-bm)/(2*(N^2-1));
  taumoptz=(-1)^N*(ap^2+2*(ap+bp)^2)/((2*(ap+bp)*(am+bm)-am*ap)*gm+(ap^2+2*(ap+bp)^2)*hm)
  taupoptz=-(am^2+2*(am+bm)^2)/((am^2+2*(am+bm)^2)*gp+(2*(am+bm)*(ap+bp)-am*ap)*hp)
  taumx=taumoptx; taupx=taupoptx; taumy=taumopty; taupy=taupopty; taumz=taumoptz; taupz=taupoptz;
  [Dx,x] = cheb(Nx); Dx=flipud(fliplr(Dx)); x=flipud(x); Dxx=Dx^2; eyex=eye(Nx+1);
  Dxx(1,:)=Dxx(1,:)-taumx*(amx*eyex(1,:)-bmx*Dx(1,:)); 
  Dxx(end,:)=Dxx(end,:)-taupx*(apx*eyex(end,:)+bpx*Dx(end,:));
  [Xeig,Sigma]=eig(Dxx.'); sigma=diag(Sigma); Xeiginv=inv(Xeig);
  [Dy,y] = cheb(Ny); Dy=flipud(fliplr(Dy)); y=flipud(y); Dyy=Dy^2; eyey=eye(Ny+1);
  Dyy(1,:)=Dyy(1,:)-taumy*(amy*eyey(1,:)-bmy*Dy(1,:)); 
  Dyy(end,:)=Dyy(end,:)-taupy*(apy*eyey(end,:)+bpy*Dy(end,:));
  [Yeig,Lambda]=eig(Dyy); lambda=diag(Lambda); Yeiginv=inv(Yeig);
  [Dz,z] = cheb(Nz); Dz=flipud(fliplr(Dz)); z=flipud(z); Dzz=Dz^2; eyez=eye(Nz+1);
  Dzz(1,:)=Dzz(1,:)-taumz*(amz*eyez(1,:)-bmz*Dz(1,:)); 
  Dzz(end,:)=Dzz(end,:)-taupz*(apz*eyez(end,:)+bpz*Dz(end,:));
  [Zeig,Omega]=eig(Dzz.'); omega=diag(Omega); Zeiginv=inv(Zeig);
  [Sigma,Lambda,Omega]=meshgrid(sigma,lambda,omega);
  coef=Sigma+Lambda+Omega;
  [X,Y,Z]=meshgrid(x,y,z); 
  ue=sin(4*pi*X).*sin(4*pi*Y).*sin(4*pi*Z);
  uex=4*pi*cos(4*pi*X).*sin(4*pi*Y).*sin(4*pi*Z);
  uey=4*pi*sin(4*pi*X).*cos(4*pi*Y).*sin(4*pi*Z);
  uez=4*pi*sin(4*pi*X).*sin(4*pi*Y).*cos(4*pi*Z);
  uexx=-16*pi^2*sin(4*pi*X).*sin(4*pi*Y).*sin(4*pi*Z);
  ueyy=-16*pi^2*sin(4*pi*X).*sin(4*pi*Y).*sin(4*pi*Z);
  uezz=-16*pi^2*sin(4*pi*X).*sin(4*pi*Y).*sin(4*pi*Z);
  lapu=uexx+ueyy+uezz; 
  lapu(:,1,:)=lapu(:,1,:)-taumx*(amx*ue(:,1,:)-bmx*uex(:,1,:)); 
  lapu(:,end,:)=lapu(:,end,:)-taupx*(apx*ue(:,end,:)+bpx*uex(:,end,:));
  lapu(1,:,:)=lapu(1,:,:)-taumy*(amy*ue(1,:,:)-bmy*uey(1,:,:)); 
  lapu(end,:,:)=lapu(end,:,:)-taupy*(apy*ue(end,:,:)+bpy*uey(end,:,:));
  lapu(:,:,1)=lapu(:,:,1)-taumz*(amz*ue(:,:,1)-bmz*uez(:,:,1)); 
  lapu(:,:,end)=lapu(:,:,end)-taupz*(apz*ue(:,:,end)+bpz*uez(:,:,end));
  tmp=permute(lapu,[3 1 2]); tmp=reshape(tmp,(Nz+1)*(Ny+1),Nx+1); tmp=tmp*Xeig; tmp=reshape(tmp,Nz+1,Ny+1,Nx+1);
  tmp=permute(tmp,[2 3 1]); tmp=reshape(tmp,(Ny+1)*(Nx+1),Nz+1); tmp=tmp*Zeig; tmp=reshape(tmp,Ny+1,(Nx+1)*(Nz+1));
  tmp=(Yeiginv*tmp); tmp=reshape(tmp,Ny+1,Nx+1,Nz+1); bb=tmp./coef;
  tmp=permute(bb,[3 1 2]); tmp=reshape(tmp,(Nz+1)*(Ny+1),Nx+1); tmp=tmp*Xeiginv; tmp=reshape(tmp,Nz+1,Ny+1,Nx+1);
  tmp=permute(tmp,[2 3 1]); tmp=reshape(tmp,(Ny+1)*(Nx+1),Nz+1); tmp=tmp*Zeiginv; tmp=reshape(tmp,Ny+1,(Nx+1)*(Nz+1));
  tmp=(Yeig*tmp); u=reshape(tmp,Ny+1,Nx+1,Nz+1); 
  err=abs(u-ue);
  maxerr=max(err(:));
  errsq=err.^2; 
  errsqocw=errsq./cwmat;
  errl2=sqrt(sum(errsqocw(:))*pi^3/(Nx*Ny*Nz));
  errl2vec=[errl2vec; errl2]; errinfvec=[errinfvec; maxerr];
  end
  format short e
  disp('L2-error   L-infinity-error'); 
  [errl2vec errinfvec]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % CHEB  compute D = differentiation matrix, x = Chebyshev grid
  function [D,x] = cheb(N)
  if N==0, D=0; x=1; return, end
  x = cos(pi*(0:N)/N)'; 
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';                  
  D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
  D  = D - diag(sum(D'));                 % diagonal entries