%  Master Program

clear;

%  Data for Model Program

xy = [ 0 0; 1 0; 0 1; 2 0; 1 1; 0 2; 2 1; 1 2; 2 2];
cn = [ 1 5 3; 1 2 5; 2 4 5; 4 7 5; 3 5 6; 5 8 6; 5 9 8; 5 7 9];
% fs = [ 4 1 500; 7 1 1000; 9 1 500];
fs = [ ];
ds = [ 1 1; 1 2; 2 2; 3 1; 4 2; 6 1];
ps = [ 9 8 1000; 8 6 1000];
E = 1e7;
nu = 0.3;
t = 1.;


%  Dimensions

nelem = size(cn,1);
nnode = size(xy,1);
nload = size(fs,1);
npres = size(ps,1);
nbc = size(ds,1);

%  Prepare the Global Stiffness Matrix

K(2*nnode,2*nnode)=0;

for ne = 1:nelem

    % Calculate the element stiffness matrix
    
    Ke = Stiffness(ne,E,nu,t,cn,xy);
    
    % Assemble into the global matrix
    
    for i = 1:3
        m = cn(ne,i);
        for j = 1:3
            n = cn(ne,j);
            K(2*m-1,2*n-1)=K(2*m-1,2*n-1)+Ke(2*i-1,2*j-1);
            K(2*m,2*n-1)=K(2*m,2*n-1)+Ke(2*i,2*j-1);
            K(2*m-1,2*n)=K(2*m-1,2*n)+Ke(2*i-1,2*j);
            K(2*m,2*n)=K(2*m,2*n)+Ke(2*i,2*j);
        end % for
    end %for
    
end %for

%  Create the load vector

F(2*nnode)=0;

%  Nodal Point forces

for nl =1:nload
    nod = fs(nl,1);
    dof = fs(nl,2);
    ndof = 2*(nod-1)+dof;
    val = fs(nl,3);
    F(ndof)=val;
end
F= F';

%  Pressures

for np = 1:npres
    n1 = ps(np,1);
    n2 = ps(np,2);
    p  = ps(np,3);
    dx = xy(n2,1)-xy(n1,1);
    dy = xy(n2,2)-xy(n1,2);
    L  = sqrt(dx^2+dy^2)
    Fx = -dy*t*p/2;
    Fy = dx*t*p/2;
    F(2*n1-1) = F(2*n1-1)+Fx;
    F(2*n1) = F(2*n1)+Fy;
    F(2*n2-1) = F(2*n2-1)+Fx;
    F(2*n2) = F(2*n2)+Fy;
end
    
    

%  Apply the displacement boundary conditions with a penalty

big = 1.e20;

for nd = 1:nbc
    nod = ds(nd,1);
    dof = ds(nd,2);
    ndof =2*(nod-1)+dof;
    K(ndof,ndof)=big;
end

% Solve for the displacements

u = K \ F;

% Get the stresses - loop over elements

for ne=1:nelem
    S(ne,:) = Stress(ne,E,nu,t,cn,xy,u);
end

