% function [p,e,tu]=tofour(p,e,t); 
%  JvP, March -97
% divide triangle(s) to four smaller ones by 
% bisecting each side. 


function [p,e,tu]=tofour(p,e,t); 

tu=[]; 

n = size(p,2); 
newp = spdiags((1:n)',0,n,n);  % vanhojen sivujen udet keskipisteet

for ti=1:size(t,2); % käy läpi kaikki kolmiot (vanhat)
   
%   jaa sivut ja  muodosta uudet pisteet 
%
%   In the point matrix P, the first and second rows contain
%   x- and y-coordinates of the points in the mesh.

  corners = t(1:3,ti);
  for ppari = [t(1:3,ti), t([2:3,1],ti)]'

% jaa edge, jos ei jaettu aikaisemmin

    if ~newp(ppari(1),ppari(2))
       n=n+1; 
       p(:,n)=p(:,ppari)*[1;1]/2; 
       newp(ppari(1),ppari(2)) = n; 
       newp(ppari(2),ppari(1)) = n; 

%  paivita (mahdollinen) muuttunut ulkosivu:
%
%        In the edge matrix E, the first and second rows contain indices
%        of the starting and ending point, the third and fourth rows 
%        contain the starting and ending parameter values, the fifth row 
%        contains the boundary segment number, and the sixth and seventh 
%        row contain the left- and right-hand side subdomain numbers.

       i1 = find(all((ppari(:,ones(1,size(e,2)))==e(1:2,:))));
       if i1~=[],
         e=e(:,[1:i1-1,i1,i1,i1+1:size(e,2)]);
         e(2,i1  )=n; 
         e(1,i1+1)=n; 
         e(3:4,i1); 
         e(3:4,i1+[0:1])=[ans(1),mean(ans);mean(ans),ans(2)];
       end
    end,

% lisätään sivu joka tapaksessa. 

    corners = [corners;newp(ppari(1),ppari(2))]; 

  end %  of if 
 
% lisää vielä uudet osakolmiot (4 kpl)

  for i = [1,4,6; 2,5,4; 4,5,6; 3,6,5].'; % suhteellinen numerointi 
    tu(:,size(tu,2)+1) = [corners(i);t(4,ti)]; 
  end

end % of for ti 
end % of funtion 
