% Script to assemble matrices for a 2D diffusion problem
% Runs in pmode - the spmd / end lines in assemb_co.m must be commented

  tic
%% set parameter values and assemble arrays
  param = p_data();
  [M1, M2, F, b,  x, e_conn] =  assemb_co(param);
  toc
  

%% Steady state solutions
  z_tmp = -full(M2)\full(F+b); % use direct solver
  z_ss  = gather(z_tmp, 1);
  toc

%% Plot and save a surface plot
  if labindex == 1
      xx = x(1:param.nodesx, 1);
      yy = x(1:param.nodesx:param.nodesx*param.nodesy, 2);
      figure
      surf(xx, yy, reshape(z_ss, param.nodesx, param.nodesy)' );
      xlabel('\bf x'); ylabel('\bf y'); zlabel('\bf T')
      t_axis = axis;
      print -dpng fig_ss.png
      close all
  end
%   
%% %% Show evolution of the Temperature
%  
  n_steps = 2; mFb = -(F+b);
  z0  = zeros(size(z_tmp), codistributor());
  M12 = full(M1 + M2);
  for jj=1:n_steps
    zn = M12\full(M1*z0 + mFb);
    z0 = zn; zg = gather(zn, 1);
    if labindex == 1     
      figure
      surf(xx, yy, reshape(zg, param.nodesx, param.nodesy)' );
      axis(t_axis); 
      print('-dpng',['fig_' int2str(jj) '.png'])
      close all
    end
  end
