% A 3D inhomogeneous Helmholtz case (u=x^2*sin(x)cos(y)cos(z)) with 
% all Dirichlet conditions
% by indirect BKM

% For different geometries, we only need change two files for boundary and test nodes.
% 1) input boundary node function, 2) test node function. 

  clear
  
  disp('Inhomogeneous 3D Helmholtz case by indirect BKM')
  
   n=input('input node number of one direction  '); %the node number of each direction
  
  t=clock;  % set initial time
  
  a=1; % length of cube
  b=1; % width of cube
  c=1; % hight of cube

   [Rx,Ry,Rz,mb]=Bcube(a,b,c,n); % boundary node generation
  
  [Dx,Dy,Dz,md]=Dcube(a,b,c,n);
  
% Lamda value of Helmholtz equation

  Lamda=sqrt(3);

% preallocating memory, mb is the total number of edge nodes
  
  Q=zeros(mb,mb); 
  bb=zeros(mb,1);   

  PQ=zeros(md,md); % preallocating memory
  db=zeros(md,1);
     
 % ---------- particular solution --------------------
 
 for i=1:md
     db(i,1)=2*sin(Dx(i))*cos(Dy(i))*cos(Dz(i))+4*Dx(i)*cos(Dx(i))*cos(Dy(i))*cos(Dz(i));
     for j=1:md
         r=D3r(Dx(i),Dx(j),Dy(i),Dy(j),Dz(i),Dz(j));  %  Euclidean distance
         PQ(i,j)=hhelm3d(r,Lamda,1); % coefficient matrix for particular solution
     end
 end

DRMNodes=md
condnumer=cond(PQ)  
alphaP=PQ\db; % evaluating the coefficients of particular solution
clear PQ db;

 % particular solutions at Dirichlet edge
 for i=1:mb
     up=0;
     for j=1:md
         r=D3r(Rx(i),Dx(j),Ry(i),Dy(j),Rz(i),Dz(j));  %  Euclidean distance
         up=up+alphaP(j)*hhelm3d(r,Lamda,2); % particular solutions at edge
     end
     bb(i,1)=Rx(i)^2*sin(Rx(i))*cos(Ry(i))*cos(Rz(i))-up;
  end

  % ----------  homogeneous solution --------------------
  
 for i=1:mb
     for j=1:mb
         r=D3r(Rx(i),Rx(j),Ry(i),Ry(j),Rz(i),Rz(j));  %  Euclidean distance
         Q(i,j)=hhelm3d(r,Lamda,0); % Dirchelet boundary entries
     end
end

mb
condnumer=cond(Q)
alpha=Q\bb; % solving the RBF equation
clear Q bb;

% evaluating some BKM solutions in boundary and domain

  [xt,yt,zt,n_test]=Dcube(a,b,c,n); % node co-ordinates of interests
  AerL2=0; % absolute error
  RerL2=0;  % relative error
 for i=1:n_test
     au=0;
     for k=1:mb  % homogeneous solutions
         r=D3r(xt(i),Rx(k),yt(i),Ry(k),zt(i),Rz(k));
         au=au+alpha(k)*hhelm3d(r,Lamda,0);
     end
     for k=1:md   % particular solutions
         r=D3r(xt(i),Dx(k),yt(i),Dy(k),zt(i),Dz(k));
         au=au+alphaP(k)*hhelm3d(r,Lamda,2); 
     end     
     u=xt(i)^2*sin(xt(i))*cos(yt(i))*cos(zt(i));  % Analytical solutions
     Aerr=u-au;  % absolute error
     if(abs(u)<10^-3)  % relative error for aboslute value <0.001
         Rerr=au;  % relative error
     else
         Rerr=(u-au)/u;  % relative error
     end
     AerL2=AerL2+Aerr^2; % square sum of absolute errors 
     RerL2=RerL2+Rerr^2; % square sum of relative errors
 end
 
AerL2=sqrt(AerL2/n_test)   % L2 absolute error
RerL2=sqrt(RerL2/n_test)   % L2 relative error

etime(clock,t)  % estimate the elapsed computing time (seconds)
InNodeNumber=(n-2)^3