% plot stress field on a regular mesh computed by ParaDiS
clc
clear all;

data1 = load('TotStress0_FMMON_4cell.out'); 
data2 = load('TotStress_pbc_FMMON_4cell.out');

ND = length(data1);
N = sqrt(ND);  M = N;
ind_shift = [floor(N/2)+1:N,1:floor(N/2)];
x = reshape(data1(:,1),N,M);
y = reshape(data1(:,2),N,M);
z = reshape(data1(:,3),N,M);
zbound=5*10^-2;

if (1)  % 2-D plot for stress along one line
    for i=1:6,
        stress1 = reshape(data1(:,i+3),N,M);
        stress2 = reshape(data2(:,i+3),N,M);
        ind_plot=floor(N/2)+0;
        figure(1); subplot(2,3,i); 
        plot([1:M],stress1(ind_plot,:),'b-',[1:M],stress2(ind_plot,ind_shift),'r--');
        %plot([1:M],stress1(ind_plot,:),'b-',[1:M],stress2(ind_plot,:), 'r--');
        xlim([1 M]);
        ylim([-zbound zbound]);
        set(gca,'FontSize',20,'FontName','Times')
        grid on
    end
    subplot(2,3,1), title('\sigma_{rr}', 'Fontsize',30);
    subplot(2,3,2), title('\sigma_{qq}', 'Fontsize',30);
    subplot(2,3,3), title('\sigma_{zz}', 'Fontsize',30);
    subplot(2,3,4), title('\sigma_{qz}', 'Fontsize',30);
    subplot(2,3,5), title('\sigma_{rz}', 'Fontsize',30);
    subplot(2,3,6), title('\sigma_{rq}', 'Fontsize',30);
end



zbound=5*10^9;

if (0) % 3-D plot for stress 
    if (1) % Print out stress1
        for i=1:6,
            stress1 = reshape(data1(:,i+3),N,M);
            figure(2); 
            subplot(2,3,i);
            mesh(stress1);
            xlim([0 M]); ylim([0 N]); zlim([-zbound zbound]);
        end
    end
    % Plot confituration for figure (2)
    subplot(2,3,1), title('(Sigma)rr');
    subplot(2,3,2), title('(Sigma)qq');
    subplot(2,3,3), title('(Sigma)zz');
    subplot(2,3,4), title('(Sigma)qz');
    subplot(2,3,5), title('(Sigma)rz');
    subplot(2,3,6), title('(Sigma)rq');
    
    
    if (1) % Print out stress2
        for j=1:6,
            stress2 = reshape(data2(:,j+3),N,M);
            figure(3);
            subplot(2,3,j);
            mesh(stress2);
            xlim([0 M]); ylim([0 N]); zlim([-zbound zbound]);
        end
    end
    % Plot confituration for figure (3)
    subplot(2,3,1), title('(Sigma)rr');
    subplot(2,3,2), title('(Sigma)qq');
    subplot(2,3,3), title('(Sigma)zz');
    subplot(2,3,4), title('(Sigma)qz');
    subplot(2,3,5), title('(Sigma)rz');
    subplot(2,3,6), title('(Sigma)rq');
end
