clear
mesh=load('mesh.dat');
segment_mesh=load('boundary_segment_mesh.dat');
u=load('u.dat');
vertices=load('vertices.dat');
stress=load('stress.dat');

s=size(vertices);
num_verts=s(1);
mesh_size=size(mesh);

node_stress=zeros(num_verts,1);

for t=1:mesh_size(1)
    s_norm=sqrt(stress(t,1)*stress(t,1)+stress(t,2)*stress(t,2)+stress(t,3)*stress(t,3)+stress(t,4)*stress(t,4));
    for i=1:3
        node_stress(mesh(t,i))=node_stress(mesh(t,i))+s_norm/3;
    end
end

for i=1:num_verts
    x(i)=vertices(i,1);
    y(i)=vertices(i,2);
    
    x_deformed(i)=x(i)+u(2*i-1);
    y_deformed(i)=y(i)+u(2*i);
end

figure
s=size(segment_mesh);
num_segments=s(1);
for s=1:num_segments
    xs(1)=x(segment_mesh(s,1));
    xs(2)=x(segment_mesh(s,2));
    
    ys(1)=y(segment_mesh(s,1));
    ys(2)=y(segment_mesh(s,2));
        
    subplot(1,2,1)
    plot(xs,ys,'k--','LineWidth',8);
    hold on
end

subplot(1,2,1)
triplot(mesh,x,y);

s=size(segment_mesh);
num_segments=s(1);
for s=1:num_segments
    xs(1)=x(segment_mesh(s,1))+u(2*segment_mesh(s,1)-1);
    xs(2)=x(segment_mesh(s,2))+u(2*segment_mesh(s,2)-1);
    
    ys(1)=y(segment_mesh(s,1))+u(2*segment_mesh(s,1));
    ys(2)=y(segment_mesh(s,2))+u(2*segment_mesh(s,2));
        
    subplot(1,2,2)
    plot(xs,ys,'k--','LineWidth',8);
    hold on
end

subplot(1,2,2)
view(2)
trisurf(mesh,x_deformed,y_deformed,node_stress,'EdgeColor','none');