clear

h = input('Hoogte van de balk (mm) : ');
L = input('Lengte van de balk (mm) : ');
nx = input('Aantal knopen over de lengte : ');
ny = input('Aantal knopen over de hoogte : ');

Emod = input('E-modulus van het materiaal (N/mm2) : ');
poiss = input('Poisson getal van het materiaal (-) : ');

apar = Emod/(1 - poiss^2);
cmat = [apar poiss*apar 0 ; poiss*apar apar 0 ; 0 0 0.5*apar*(1-poiss)];

nne = 3;
ndf = 2;
q = 1;          % Belasting in N/mm

% Read finite element mesh
[nnm,nem,xg,con] = makemesh(nx,ny,h,L);

% Read boundary conditions
[ndisp,nload,dispval,dispfix,forcval,forcfix] = dispvalues(nnm,nem,con,h,L,xg,q);

% Prepare the global stiffness matrix
[jdiag,ntrm] = prepare(nem,nnm,nne,ndf,con);

kglob = zeros(ntrm,1);
fglob = zeros(nnm*ndf,1);
    

for ielm = 1:nem
    xloc(:,:) = xg(con(ielm,:),:);
    [nu,bu,detJ] = Tri3(xloc,nne,ndf);
    kmat = bu'*cmat*bu*detJ;
    [kglob] = assemble(nne,ndf,con(ielm,:),jdiag,kmat,kglob); 
end


for i=1:ndisp
    [kglob,fglob] = boundary(nnm*ndf,jdiag,dispval(i,1),dispfix(i,1),kglob,fglob);
end

for i=1:nload
    fglob(forcfix(i,1),1) = fglob(forcfix(i,1),1) + forcval(i,1);
end

[kglob] = reduce(nnm*ndf,ntrm,jdiag,kglob);
[fglob] = solve_direct(nnm*ndf,ntrm,jdiag,kglob,fglob);

for i=1:nnm
    usol(i,1) = fglob(2*i-1,1);
    usol(i,2) = fglob(2*i,1);
end

% Compute reaction forces and nodal values of strains and stresses

strain = zeros(nnm,3);
stress = zeros(nnm,3);
area = zeros(nnm,1);

for ielm = 1:nem
    
    xloc(:,:) = xg(con(ielm,:),:);
    [nu,bu,detJ] = Tri3(xloc,nne,ndf);
    
    for i=1:nne
        uloc(2*i-1,1) = usol(con(ielm,i),1);
        uloc(2*i,1)   = usol(con(ielm,i),2);
    end

    strainloc = bu*uloc;
    stressloc = cmat*bu*uloc;
    for i=1:nne
        for j=1:3
            strain(con(ielm,i),j) = strain(con(ielm,i),j) + strainloc(j,1)*detJ/3;
            stress(con(ielm,i),j) = stress(con(ielm,i),j) + stressloc(j,1)*detJ/3;
        end
        area(con(ielm,i),1) = area(con(ielm,i),1) + detJ/3;
    end
       
end

for i=1:nnm
    stress(i,:) = stress(i,:)/area(i,1);
    strain(i,:) = strain(i,:)/area(i,1);
end

% trisurf(con,xg(:,1),xg(:,2),stress(:,1));
% view(0,90);
% axis equal;