function [pnew,tnew ] = mesh2acf( pinput,tinput )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MESH2ACF
%% the function i) converts [p,t] produced by DISTMESH2D  
%% to the format required by ACF software, ii) detects boundary, 
%% and iii) refines the grid if desired
%% It makes small changes to the grid, e.g, unifies orientation
%% It prints <coordinates.dat>, <elemenst3.dat>, <boundary.dat>
%% It can also execute (uniform) grid refinement 
%% in which case it outputs pnew and tnew
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% INPUT: pinput,tinput: points and triangles from any triangulation
%%% OUTPUT: prints files <coordinates.dat>, <elemenst3.dat>, <boundary.dat>
%%% if output requested, the code refines grud iniformly to <pnew,tnew>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% M. Peszynska, 2013
%% © Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if nargout > 1, Ifrefine =1; else Ifrefine =0; end
    
    poutput =  zeros(size(pinput,1),size(pinput,2)+1);
    toutput =  zeros(size(tinput,1),size(tinput,2)+1);
    toutput(:,2:4) = tinput(:,1:3);
    %%%% 
    MAXTRISPERNODE = 2;
    Edges = zeros(size(tinput,1)*3,2 + MAXTRISPERNODE);
    Lastedge = 0;

    %%% records edge numbers for each triangle
    %% Triangle_edge(j,1) is edge 1 of Triangle j 
    Triangle_edge = 0*tinput;
    
     function add_edge(in1,in2,j)
        %% fprintf('Add edge [%d %d] with tri %d\n',in1,in2,j);
        n1=min(in1,in2);n2=max(in1,in2); 
        record=0;indl=[];indr=[];
         if Lastedge == 0, 
            record=1;
         else
            ledge = Edges(1:Lastedge,1);
            redge = Edges(1:Lastedge,2);
            indl=find(ledge==n1);
            if length(indl)==0,record=1;
            else
                indr=find(redge(indl)==n2);
                if length(indr)==0,record=1;end;
            end
         end
         %% fprintf('....indl=%d indr =%d record=%d\n',length(indl),length(indr),record);
        if record==1
            Lastedge = Lastedge+1;
            Edges(Lastedge,1)=n1;Edges(Lastedge,2)=n2;
            Edges(Lastedge,3) = 1;
            Edges(Lastedge,4)=j;
        else
           if length(indr) >1, error('too many edges'); end;
           indexedge = indl(indr);
           %%% current number of triangles to which an edge belongs
           ntris = Edges(indexedge,3);
           % fprintf('...found %d old tris for edge=%d\n',ntris,indexedge);
           if ntris<MAXTRISPERNODE
                Edges(indexedge,3) = ntris+1;
                Edges(indexedge,3+ntris+1)=j;
           else %% grow matrix and record
               error('must grow matrix');
           end
        end
     end
    
    %%%%%%%%%%%%%%%%
    minx=min(pinput(:,1));maxx=max(pinput(:,1));
    miny=min(pinput(:,2));maxy=max(pinput(:,2));
    
    delx = maxx-minx; dely = maxy-miny;
    ave_area = dely*delx/size(tinput,1);
    SCALE = 0.1;
    bbox = [minx-delx*SCALE, maxx+delx*SCALE, miny-dely*SCALE, maxy+dely*SCALE]; 
    
    clf;hold on;
    for j=1:size(pinput,1)
        poutput(j,2:3) = pinput(j,1:2);
        poutput(j,1) = j;
        text(pinput(j,1),pinput(j,2),sprintf('%d',j));
        axis(bbox);
    end
    title('Nodes of FE grid');
    pause,
    save('coordinates.dat','poutput','-ASCII');
    
    
    for j=1:size(tinput,1)
        toutput(j,1)=j;
        tri = tinput(j,:);
        xcoor = pinput(tri,1); ycoor = pinput(tri,2);
       
        
        plot(xcoor(1:2),ycoor(1:2),'-');
        plot(xcoor(2:3),ycoor(2:3),'-');
        plot([xcoor(3) xcoor(1)],[ycoor(3) ycoor(1)],'-');
        text(sum(xcoor)/3,sum(ycoor)/3,sprintf('%d',j));
        title(sprintf('Triangle %d <%d %d %d>',j,tri));
        pause(0.05);
        
        areatri = [xcoor, ycoor, ones(3,1)];
        dtri = det(areatri);
        %% fprintf('Triangle %d <%d %d %d> area=%g\n',j,tri,1/2*dtri);
        if abs(dtri)<1e-3*ave_area,error('degenerate triangle');end

        if dtri < 0 %% reverse orientation
            toutput(j,3)=tri(3);
            toutput(j,4)=tri(2);
            fprintf('triangle %d fixed. From %d %d %d to %d %d %d \n',j,...
                tri,toutput(j,2:4));
        end
        %%% create a list of edges to identify the boundary edges
        if 1
        add_edge(toutput(j,2),toutput(j,3),j);
        add_edge(toutput(j,3),toutput(j,4),j);
        add_edge(toutput(j,4),toutput(j,2),j);
        end    
        %%% find which edge has only one triangle 
    end
    save('elements3.dat','toutput','-ASCII');
    
   %%%%%%% show the edges 
   fprintf('Number of edges=%d\n',Lastedge);
   %% Edges(1:Lastedge,:),pause,
               
   %%% show for each triangle the edges that form it 
   %% this is only necessary for quick lookup and mesh refinements
   for j=1:Lastedge
       %%% loop over each of the (possibly) two triangles for this edge
       for k=1:Edges(j,3)  %% one or two
           %% Edges(j,4) is first triangle number and Edges(j,5) is second traingle
           mytrinumber = Edges(j,3+k);
           if Triangle_edge(mytrinumber,1)==0
               Triangle_edge(mytrinumber,1)=j;
           elseif Triangle_edge(mytrinumber,2)==0
               Triangle_edge(mytrinumber,2)=j;
           elseif Triangle_edge(mytrinumber,3)==0
               Triangle_edge(mytrinumber,3)=j;
           else
               error('too many edges ?');
           end
       end
   end
   %% fprintf('Edges for each triangle');
   %% Triangle_edge,
   %% pause,
   
   %%% identify and print boundary edges
   indb = find(Edges(1:Lastedge,3)==1);
   boundary=zeros(length(indb),3);
   for j=1:length(indb)
       boundary(j,1)=j;
       boundary(j,2:3)=Edges(indb(j),1:2);
   end
   save('boundary.dat','boundary','-ASCII');
   hold off;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uniform grid refinement
   % Ifrefine = 1;
   
   if Ifrefine ==0, return; end
   
   oldnodes = size(pinput,1);
   pnew = zeros(size(pinput,1)+Lastedge,2);
   pnew(1:size(pinput,1),1:2) = pinput(1:end,1:2);
   
   %%% new nodes added at midpoint of each edge
   for j=1:Lastedge
       n1=Edges(j,1);n2=Edges(j,2);
       pnew(oldnodes+j,1) = 0.5*(pinput(n1,1)+pinput(n2,1));
       pnew(oldnodes+j,2) = 0.5*(pinput(n1,2)+pinput(n2,2));
   end
   %pnew,pause 
   %%% new triangles are formed as quarters of the original triangle

   tnew = zeros(size(tinput,1)*4,3);
   k = 0;
   for j=1:size(tinput,1)
       %%% old nodes of the triangle j
       mynodes = tinput(j,1:3);
       
       myedges = Triangle_edge(j,1:3);
       me12 = myedges(1); %%% me23 = mydeges(2); me31 = myedges(3); 
       
       %% identify edges with nodes
       me1 = Edges(me12,1); me2 = Edges(me12,2); me3 = setdiff(mynodes,[me1,me2]);
       
       %% find me13 between nodes ne1 and ne3
       if Edges(myedges(2),1) == me1 || Edges(myedges(2),2) == me1
           me13 = myedges(2);
       else
           me13 = myedges(3);
       end
       me23 = setdiff(myedges,[me12, me13]);
       
       %% the midpoint nodes actually have numbers 
       me12 = me12 + oldnodes; me23 = me23 + oldnodes; me13 = me13 + oldnodes;
       
       %% now we have identified 6 nodes (3 old 3 new) in triangle j
       %% create new trainagles
       k = k+1;
       tnew (k,1) = me1; tnew(k,2) = me12; tnew(k,3)=me13;
       
       k = k+1;
       tnew (k,1) = me12; tnew(k,2) = me2; tnew(k,3)=me23;
       
        k = k+1;
       tnew (k,1) = me23; tnew(k,2) = me3; tnew(k,3)=me13;
       
        k = k+1;
       tnew (k,1) = me12; tnew(k,2) = me23; tnew(k,3)=me13;
       
       
   end
   
end         
