% -------------------------------------------------------------------------------------- %
% Methodes de discretisation en fluides                                                  %
% M.A. Habisreutinger                                                                    %
% Vendredi 10 avril 2019                                                                 %
% -------------------------------------------------------------------------------------- %

function SteadyPoiseulle3D()
  % M.A. Habisreutinger 2019
  
  %% 1.   PARAMETERS

  %  1.1  Mesh
  ndimensions = 3 ;
  L = [8, 1, 2] ;
  npoints = [8, 16, 8] ;

  %  1.2  Diffusion coefficient
  mu = 1e+00 ;
  
  %  1.3  Pressure drop
  dpdx = -1e+00 ;

  %  1.4  Vtk output for Paraview (download at www.paraview.org)
  write = true ;
  filename = 'Poiseuille3D' ;


  %% 2.   COMPUTATION
  
  %  2.1  Number of variables, nodes and unknowns
  nvars = ndimensions+1 ;
  nnodes = prod(npoints) ;
  nuks = nvars*nnodes ;

  %  2.1  Nodal distributions 
  for idim = 1 : ndimensions
    xi{idim} = linspace(-L(idim)/2, +L(idim)/2, npoints(idim)) ;
    h(idim) = L(idim)/(npoints(idim)-1) ;
  end

  %  2.2  Mesh nodes
  [X{1}, X{2}, X{3}] = ndgrid(xi{1}, xi{2}, xi{3}) ;
  x = zeros(nnodes, ndimensions) ;
  for idim = 1 : ndimensions 
    x(:, idim) = reshape(X{idim}, nnodes, 1) ;
  end
  
  %  2.4  Boundary descriptors
  left = x(:, 1) == -L(1)/2 ;
  right = x(:, 1) == +L(1)/2 ;
  bottom = x(:, 2) == -L(2)/2 ;
  top = x(:, 2) == +L(2)/2 ;
  back = x(:, 3) == -L(3)/2 ;
  front = x(:, 3) == +L(3)/2 ;

  %  2.3  Discrete operators
  for idim = 1 : ndimensions
    e = ones(npoints(idim), 1) ;
    % Central finite differences (2nd order)
    Mi{idim} = spdiags(e, 0, npoints(idim), npoints(idim)) ;
    Ai{idim} = spdiags(-mu*(1/h(idim)^2)*[+1*e -2*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
    Di{idim} = spdiags(+ (1/(2*h(idim)))*[-1*e  0*e +1*e], [-1, 0, +1], npoints(idim), npoints(idim)) ;
    % Neumann boundary conditions
    Ai{idim}(1, 1) = -mu*(-1/h(idim)^2) ;
    Ai{idim}(end, end) = -mu*(-1/h(idim)^2) ;
    % Shifted scheme (1st order)
    Di{idim}(1, 1) = -1/h(idim) ;
    Di{idim}(1, 2) = +1/h(idim) ;
    Di{idim}(end, end-1) = -1/h(idim) ;
    Di{idim}(end, end) = +1/h(idim) ;
  end

  %  2.4  Tensorisation for a 3D problem
  Mo = kron(Mi{3}, kron(Mi{2}, Mi{1})) ;
  Ao = kron(Ai{3}, kron(Mi{2}, Mi{1})) ...
     + kron(Mi{3}, kron(Ai{2}, Mi{1})) ...
     + kron(Mi{3}, kron(Mi{2}, Ai{1})) ;
  Dx = kron(Mi{3}, kron(Mi{2}, Di{1})) ;
  Dy = kron(Mi{3}, kron(Di{2}, Mi{1})) ;
  Dz = kron(Di{3}, kron(Mi{2}, Mi{1})) ;
  O = sparse(nnodes, nnodes) ;

  %  2.7  Assembling global matrices
  M = [Mo, O,  O, O; ...
        O, Mo, O, O; ...
        O, O, Mo, O; ...
        O, O,  O, O] ;

  A = [Ao, O,  O,  Dx; ...
        O, Ao, O,  Dy; ...
        O, O,  Ao, Dz; ... 
       Dx, Dy, Dz, O] ;
 
  %  2.6  Forcing term
  f = zeros(nuks, 1) ;
  Mf = M*f ;
  
  %  2.5  Exact solution
  vx = (1/2)*dpdx*(x(:, 2)-L(2)/2).*(x(:, 2)+L(2)/2) ;
  vy = zeros(nnodes, 1) ;
  vz = zeros(nnodes, 1) ;
  p = dpdx*x(:, 1) ;
  u = [vx; vy; vz; p] ;
  
  %  2.10 Correction for Dirichlet boundary conditions
  c = [(bottom | top | back | front); ...
    +  (bottom | top | left | right | back | front); ...
    +  (bottom | top | left | right | back | front); ...
    +  (left | right)] ;

  C = spdiags( c, 0, nuks, nuks) ;
  F = spdiags(~c, 0, nuks, nuks) ;

  Mf = F*Mf+C*u ;
  A = F*A+C ;
 
  %  2.13 Solution of the discrete equations (direct method)
  uh = A\Mf ;

  %  2.14 Retrieving the variables
  indices = 1:nnodes ;
  vxh = uh(indices) ;
  vyh = uh(indices+nnodes) ;
  vzh = uh(indices+2*nnodes) ;
  ph = uh(indices+3*nnodes) ;


  %% 3.   DATA ANALYSIS

  %  3.1  Error norms
  eh = norm(M*(uh-u))/norm(M*u) ;

  %  3.2  Display
  disp(['Error norm : ', num2str(eh)]) ;

  %  3.3  Global matrices plot
  figure(1) ;
  spy(A) ;
  xlabel('j') ; ylabel('i') ; 
  title('Structure of the global discretization matrix') ;
  
  %  3.2  Vtk output (visualisation tool kit à lire dans paraview)
  if write
    WriteVtkOutput(filename, {'vx', 'vy', 'p'}, xi, {vxh, vyh, ph}, npoints) ;
  end

end

% -------------------------------------------------------------------------------------- %

function WriteVtkOutput(filename, varnames, xi, fields, npoints)
  % M.A. Habisreutinger 2019
  
  % File
  file = fopen(strcat(filename, '.vtk'), 'w') ;

  % Header
  fprintf(file, '%s\n', '# vtk DataFile Version 3.1') ;
  fprintf(file, '%s\n', 'Export Q') ;
  fprintf(file, '%s\n', 'ASCII') ;
  fprintf(file, '%s\n', 'DATASET RECTILINEAR_GRID') ;
  
  fprintf(file, '%s ', 'DIMENSIONS') ;
  fprintf(file, '%u %u %u\n', npoints) ;
  
  fprintf(file, '%s ', 'X_COORDINATES') ;
  fprintf(file, '%u ', npoints(1)) ;
  fprintf(file, '%s\n', 'FLOAT') ;
  fprintf(file, '%15.5f\n', xi{1}) ;
  
  fprintf(file, '%s ', 'Y_COORDINATES') ;
  fprintf(file, '%u ', npoints(2)) ;
  fprintf(file, '%s\n', 'FLOAT') ;
  fprintf(file, '%15.5f\n', xi{2}) ;
    
  fprintf(file, '%s ', 'Z_COORDINATES') ;
  fprintf(file, '%u ', npoints(3)) ;
  fprintf(file, '%s\n', 'FLOAT') ;
  fprintf(file, '%15.5f\n', xi{3}) ;
  
  fprintf(file, '%s ', 'POINT_DATA') ;
  fprintf(file, '%u\n', prod(npoints)) ;
  
  % Fields
  field = reshape(fields{1}, prod(npoints), 1) ;
  fprintf(file, '%s ', 'SCALARS') ;
  fprintf(file, '%s ', varnames{1}) ;
  fprintf(file, '%s\n', 'FLOAT') ;
  fprintf(file, '%s\n', 'LOOKUP_TABLE default') ;
  fprintf(file, '%15.5f\n', field') ;
  
  field = reshape(fields{2}, prod(npoints), 1) ;
  fprintf(file, '%s ', 'SCALARS') ;
  fprintf(file, '%s ', varnames{2}) ;
  fprintf(file, '%s\n', 'FLOAT') ;
  fprintf(file, '%s\n', 'LOOKUP_TABLE default') ;
  fprintf(file, '%15.5f\n', field') ;
  
  field = reshape(fields{3}, prod(npoints), 1) ;
  fprintf(file, '%s ', 'SCALARS') ;
  fprintf(file, '%s ', varnames{3}) ;
  fprintf(file, '%s\n', 'FLOAT') ;
  fprintf(file, '%s\n', 'LOOKUP_TABLE default') ;
  fprintf(file, '%15.5f\n', field') ;
  
  fclose(file) ;
  
end

% -------------------------------------------------------------------------------------- %