function [t,subdomains]=domain_1d(u0,num_subdomains,size_of_subdomain,num_iter)
% DOMAIN_1D
%
% [t,subdomains]=domain_1d(u0,num_subdomains,size_of_subdomain,num_iter)
%
% Solves Laplace's Equation in one-dimension using domain
% decomposition.
%
% u0 = sets initial boundary conditions and implicitly determines the grid
% size.  u0(1) and u0(biggest_index) are the left and right Dirichlet
% boundary conditions.  the rest of u0 is ignored.  the size of the
% grid is the size of u0.
%
% the other parameters should be self-explanatory; if not, examine
% the code and comments within.
%
% RETURNS:
%
% t = the range of coordinates of each subdomain ("the domain of
%   each domain")
%
% subdomains = the Laplace solution on each subdomain.  It is up to
%   you to glue these overlapping subdomains into one complete
%   solution over the whole domain
%
% $Id: domain_1d.m,v 1.4 2004/03/31 07:48:54 kenta Exp $


% In general (in 1D), when specifying HOW to decompose a domain, there
% are four variables: total number of grid points, number of
% subdomains, the size of each subdomain, and the amount of overlap.
% Choosing any three of the four determines the fourth. This example
% code will make overlap the dependent variable; however, you are
% welcome go differently.

num_grid_points=size(u0,2);
left_edge_of_rightmost_domain=num_grid_points-size_of_subdomain;

% This won't do the right thing for num_subdomains==1
stepping= left_edge_of_rightmost_domain/(num_subdomains-1);

% Omit semicolon and output the size of the overlap.  The amount of
% overlap doesn't always come out to an integer -- in which case it is
% the average overlap.  There won't be convergence if overlap==1
% anywhere (think about it).
overlap = size_of_subdomain - stepping 


% preallocate memory
% In serial 1D code, all the subdomains fit into a nice neat
% matrix.  In 2D, and in parallel, this will be distributed data
% structure across processes.

% We zero out the internal grid points.  However, one really ought
% "precondition" the internal points with some approximate solution
% in order to speed convergence.

subdomains=zeros(num_subdomains, size_of_subdomain);

% Set the global left and right boundary conditions
subdomains(1,1)=u0(1);
subdomains(num_subdomains,size_of_subdomain)=u0(num_grid_points);


% Calculate the domain of each subdomain.
t=[];
for dom=1:num_subdomains

  left=(dom-1)*stepping+1;
  left_rounded=round(left);
  right=left_rounded+size_of_subdomain-1; 
  % the last "right" should be num_subdomains
  % if I did this right.

  t=[t;left_rounded:right];
end

% We will need the offset of each subdomain with respect to its
% neighbors for communication between domains later.
offsets=diff(t(:,1));

% Instead of hard-coding the number of iterations as is done in this
% loop, you will want to implement a test for convergence.
for iter=1:num_iter
  
  % This loop is the messy communications between processors you
  % will have to implement  
  for dom=1:num_subdomains-1
    
    % Update the right boundary condition by pulling an internal
    % value of the subdomain to the right.
    subdomains(dom,size_of_subdomain)= ...
        subdomains(dom+1,size_of_subdomain-offsets(dom));
    
    % Similarly for left
    subdomains(dom+1,1)=subdomains(dom,(1+offsets(dom)));
  end

  %temp=subdomains  % uncomment this if you want to see the
                   % intermediate update of boundary conds.
  
  % You should parallelize this loop by having each subdomain
  % solved on a different process.
  for dom=1:num_subdomains
    subdomains(dom,:)=laplace_1d(subdomains(dom,:));
  end

end

% Finally, there is the problem of gluing the overlapping subdomains
% together into one.  This is left up to you.
