function [bc,df,edof,dof,coord,enod,plot_dof] = TopConstMod_Assignment2015(p,t,dtau_x,du,th,control) 
% [bc,df,edof,dof,coord,enod,plot_dof] = TopConstMod_Assignment2015(p,t,dtau_x,du,th,control) 
%-------------------------------------------------------------
% PURPOSE
% Extract bc,df,edof,dof,coord,enod,plot_dof from a mesh generated in
% pdetool. Note that this routine is specifically written for the geometry
% and load case given in the computer lab in constitutive modelling 2014. 
%
%
% INPUT:  p         Variable retrieved from pdetool contains 
%                   x- and y- coordinates of all nodes
%         t         Variable retrieved from pdetool contains information 
%                   of what nodes are connected to each element
%         dtau_x    Incremental traction stress applied at the left boundary
%         du        Incremental displacement on the left boundary
%         th        thickness
%         control   control=0 force controlled displacement
%                   control=1 displacement controlled displacement
%
% OUTPUT: bc        First column dof, second column known displacement
%         df        External applied force each time step (ndof X 1)
%         edof      Element topology matrix (how the dofs are conected to each element) (nelm X 7)
%         dof       Dof conexted to each node (nnod X 2)
%         coord     x- and y- coordinates of all nodes (nnod X 2)
%         enod      Element topology matrix (how the nodes are conected to each element) (nelm X 3)     
%         plot_dof  Degree of freedome used for plottin the
%                   force-displacement curves (the x-dir of the lower left node)
%-------------------------------------------------------------
% LAST MODIFIED: H Askfelt 2015-10-02
%-------------------------------------------------------------

nen=3;                                  % number of nodes per element
enod=t(1:nen,:)';                       % 
nelm=size(enod,1);                      % number of elements  
coord=p';                               %
dofnode=2;                              % number of dofs per node
nnod=size(coord,1);                     % number of nodes
dof=[(1:nnod)',(nnod+1:2*nnod)'];       % give each dof a number
ndof=max(max(dof));                     % number of dofs  
edof=zeros(nelm,nen*dofnode+1);         % allocate space for edof

% Generate edof from enod and dof      
for ie=1:nelm
   edof(ie,:)=[ie dof(enod(ie,1),:), dof(enod(ie,2),:), dof(enod(ie,3),:)];
end

[ex,ey]=coordxtr(edof,coord,dof,nen);
eldraw2(ex,ey,[1 2 0])

% Find the coordinates for the corners
minx=min(coord(:,1));                   % minimum x-coord. of the geometry
maxx=max(coord(:,1));                   % maximum x-coord. of the geometry
miny=min(coord(:,2));                   % minimum y-coord. of the geometry
maxy=max(coord(:,2));                   % maximum x-coord. of the geometry

% Define bc and df
bc=[];                                  % create an empty matrix bc
df=zeros(ndof,1);                       % allocate space for df

node_minx=[];                           % allocate space for node_minx 
                                        % node_minx will contain all nodes on the left boundary
for inod=1:nnod %loop through all nodes
    
    if (abs(coord(inod,1)-minx)<1e-12)            % if the current node inod is on the left boundary
        if(control==1 || control==2 )   % if we are solving a displacement controlled problem
            bc=[bc;dof(inod,1), du];    % set the displacement in x-dir of inod to du
        end
        node_minx=[node_minx;inod];     % add inod to node_minx
        if (abs(coord(inod,2)-maxy)<1e-12)
            plot_dof=dof(inod,1);       % saving the dof in x-dir of the node in the upper left corner
        end

    elseif(abs(coord(inod,1)-maxx)<1e-12)         % if the current node inod is on the right boundary
%         inod
%         coord(inod,2)
        bc=[bc;dof(inod,1), 0];         % symmetry on right boundary no dissplacement
    end

    if(abs(coord(inod,2)-miny)<1e-12)             % if the current node inod is on the lower boundary
        bc=[bc;dof(inod,2), 0];         % symmetry on lower boundary no dissplacement
    end
end



tmp=[node_minx coord(node_minx,2)];     % create a matrix tmp containing all 
                                        % nodes on the left boundary and
                                        % there y-coordinate
[Y,I]=sort(tmp(:,2));                   % sort tmp by the nodes y-coordinates
tmpsort=tmp(I,:);
 for inod=1:size(tmp,1)
    tnod=tmpsort(inod,1);
    if inod==1  % smalest y-val on the left bc
        L=(tmpsort(2,2)-tmpsort(1,2));  % external boundary length of the element
        df(dof(tnod,1))=dtau_x*L/2;        
    elseif inod==size(tmp,1) % largest y-val on the left bc
        L=(tmpsort(size(tmp,1),2)-tmpsort(size(tmp,1)-1,2)); % external boundary length of the element
        df(dof(tnod,1))=dtau_x*L/2;
    else    %Gets contribution from both the element above and below
    %element below
    L=(tmpsort(inod,2)-tmpsort(inod-1,2));  % external boundary length of the element
    df(dof(tnod,1))=df(dof(tnod,1))+L/2*dtau_x;
    %element above
    L=(tmpsort(inod+1,2)-tmpsort(inod,2));  % external boundary length of the element
    df(dof(tnod,1))=df(dof(tnod,1))+L/2*dtau_x;
    end
 end
 
 df = df*th;                        % multiply df with the thickness
 
 

% % Below are optional plot routines where the added bc and df may be seen
% a=zeros(ndof,1);
% a(bc(:,1))=1;
% ed=extract(edof,a);
% % ed(1,:)=ed(1,:)*1e5
% figure(1)
% fill(ex',ey',ed(:,[1 3 5])');     %(markes all dofs with bc in x-dir)
% figure(2)
% fill(ex',ey',ed(:,[2 4 6])');     %(markes all dofs with bc in y-dir)
% 
% %  keyboard
% 
% a=zeros(ndof,1);
% a(bc(:,1))=bc(:,2);
% ed=extract(edof,a);
% figure(3)
% fill(ex',ey',ed(:,[1 3 5])');   %(markes all dofs with bc in x-dir)
% figure(4)
% fill(ex',ey',ed(:,[2 4 6])');   %(markes all dofs with bc in y-dir)
% 
% 
% 
% ed=extract(edof,df);
% figure(5)
% fill(ex',ey',ed(:,[1 3 5])');   %(markes all dofs with fl in x-dir)
% figure(6)
% fill(ex',ey',ed(:,[2 4 6])');   %(markes all dofs with fl in y-dir)
% 
% % keyboard



% for ie=1:nelm
%     for idof=1:3
%         
%         if (enod(ie,idof)==13)
%             ie
%             enod(ie,:)
%         end
%     
%     end
% end    
% keyboard
