
%program to plot FEM results
%function plot_FEM_results(coreDir, pixelsize,filename3, filename4)
%make sure cell flders added in FEM and Tractions
%this file writes averaged intercellular stress values as a .dat file and
%saves entire stress matrix for each time frame

%miku, jan 3,2017


function plot_FEM_results(File,Settings)

for j=1:File.NFiles.Fluorescence,
    copyfile(['C:\code\Trac_MSM_Code\Data\Tractions\Time',num2str(j)],['C:\code\Trac_MSM_Code\Data\FEM\Time',num2str(j)]);
end
    
cd ('C:\code\Trac_MSM_Code\Data');


for k=1:File.NFiles.Fluorescence,
 

% plot data
pixelTometer = Settings.PixelSize*(1e-6); 
tol = 1.0e-12; % a number smaller compared to the grid size

% read tractions
%tractions = textread([File.pathname,'Tractions','/', 'Time', num2str(k), '/', 'traction', '.dat']);
tractions = textread(['C:\code\Trac_MSM_Code\Data\','Tractions','\', 'Time', num2str(k), '\', 'traction', '.dat']);



% read node_coords
%node = textread([File.pathname, 'FEM','/', 'Time', num2str(k), '/', 'node_coords', '.dat']);
node = textread(['C:\code\Trac_MSM_Code\Data\', 'FEM','\', 'Time', num2str(k), '\', 'node_coords', '.dat']);
node = node/pixelTometer; % rescale the nodal coordinates to match with traction grid


%read results
%stresses = textread([File.pathname, 'FEM','/', 'Time', num2str(k), '/', 'stresses', '.dat']);
stresses = textread(['C:\code\Trac_MSM_Code\Data\','FEM','\', 'Time', num2str(k), '\', 'stresses', '.dat']);
plotVar = stresses(:,1:7);%includes StrsX, StrsY, StrsXY, Pmax, Pmin, x or stress, y or stress


% grid the stresses
plotVar_onTracGrid = zeros(size(tractions,1),7);

%solve for various values below

avg_norm_stress = mean((stresses(:,4)+stresses(:,5))/2) 
max_shear_stress = mean((stresses(:,4)-stresses(:,5))/2)


if (strcmpi(Settings.EvenStresses,'Yes')) == 1 



for i=1:size(node,1)
    for j=1:size(tractions,1)
        if((tractions(j,1)<node(i,1)+tol) && (tractions(j,1)>node(i,1)-tol))
            if((tractions(j,2)<node(i,2)+tol) && (tractions(j,2)>node(i,2)-tol))
                plotVar_onTracGrid(j,:) = plotVar(i,:);
                break;
            end
        end
    end
end

NewStresses.Stresses{k} = plotVar_onTracGrid;%save stresses

if (strcmpi(Settings.PlotStresses,'Yes')) == 1 %plot normal and shear stress


avg_norm = (plotVar_onTracGrid(:,4) + plotVar_onTracGrid(:,5))/2;
max_shear = (plotVar_onTracGrid(:,4) - plotVar_onTracGrid(:,5))/2;
node_x = reshape(tractions(:,1), Settings.domSize(2),Settings.domSize(1));
node_y = reshape(tractions(:,2), Settings.domSize(2),Settings.domSize(1));
stress_norm = reshape(avg_norm, Settings.domSize(2),Settings.domSize(1));%avg norm stress
stress_shear = reshape(max_shear, Settings.domSize(2),Settings.domSize(1));%max shear


        if (strcmpi(Settings.reshapeMatrix,'Yes')) == 1 
        
            %resize matrix to reduce pixelation of output
            %TracX = imresize(TracX, size(TracX)*2,'bicubic');
            %TracY = imresize(TracY, size(TracY)*2,'bicubic');   
            stress_norm = imresize(stress_norm, size(stress_norm)*2,'bicubic');
            stress_shear = imresize(stress_shear, size(stress_shear)*2,'bicubic');
            node_x=imresize(node_x, size(node_x)*2,'bicubic');
            node_y=imresize(node_y, size(node_y)*2,'bicubic');  
        
        end
   
%plot avg normal stress 
figure(1); 

   surface(node_x,node_y,stress_norm, 'edgecolor','none'); 
   view(2); colormap jet; shading flat;
   cbh = colorbar; set(cbh,'FontSize',18,'FontWeight','bold');
   set(gca,'clim',[Settings.avg_norm_min Settings.avg_norm_max])
   axis equal; axis tight; axis ij; axis off    
   axis([min(min(node_x)) max(max(node_x)) min(min(node_y)) max(max(node_y))]) 
   set(gca,'FontSize',18,'FontWeight','bold'); box on;
   title(' Average Normal Stress (Pa)','FontSize',18,'FontWeight','bold');
   
   
    if (strcmpi(Settings.height,'infinite')) == 1
          saveas(gcf,[File.pathname, 'Results', '/AvgNormStress', num2str(k), '.jpg']);
    end
    
    if (strcmpi(Settings.height,'finite')) == 1
        saveas(gcf,[File.pathname, 'Results', '/AvgNormStress_finite', num2str(k), '.jpg']);
    end
    
    if (strcmpi(Settings.height,'constrained')) == 1
        saveas(gcf,[File.pathname, 'Results', '/AvgNormStress_constrained', num2str(k), '.jpg']);
    end
    
 %plot maximum shear intercellular stress
  figure(4); 
   surf(node_x,node_y,stress_shear, 'edgecolor','none'); 
   view(2); colormap jet; shading flat; 
   cbh = colorbar; set(cbh,'FontSize',18,'FontWeight','bold');
   set(gca,'clim',[0 Settings.max_shear_max])
   axis equal; axis tight; axis ij; axis off  
   axis([min(min(node_x)) max(max(node_x)) min(min(node_y)) max(max(node_y))]) 
   set(gca,'FontSize',18,'FontWeight','bold'); box on;
   title('Maximum Shear Stress (Pa)','FontSize',18,'FontWeight','bold'); 
   
  if (strcmpi(Settings.height,'infinite')) == 1
          saveas(gcf,[File.pathname, 'Results', '/MaxShearStress', num2str(k), '.jpg']);
    end
    
    if (strcmpi(Settings.height,'finite')) == 1
        saveas(gcf,[File.pathname, 'Results', '/MaxShearStress_finite', num2str(k), '.jpg']);
    end
    
    if (strcmpi(Settings.height,'constrained')) == 1
        saveas(gcf,[File.pathname, 'Results', '/MaxShearStress_constrained', num2str(k), '.jpg']);
    end
 
    
end
    

end
   %write intercellular stress values into file
   
  
         results = [avg_norm_stress max_shear_stress];
         
         if k == 1,
             
             finalResults = results;
             
         else
             
             finalResults = [finalResults;results]
        
         end
         
         
         %save averaged stress values into .dat file
         if (k == File.NFiles.Fluorescence)%File.NFiles.Fluorescence
             
             
             if (strcmpi(Settings.height,'infinite')) == 1
         
                %f3 = ([File.pathname, '/Results','/', 'StressResults_unconstrained', '.dat']);
                %dlmwrite(f3,finalResults,'precision',6);
                f3 = (['C:\code\Trac_MSM_Code\Data\', 'Results','\', 'StressResults_unconstrained', '.dat']);
                dlmwrite(f3,finalResults,'precision',6);
                
             elseif (strcmpi(Settings.height,'finite')) == 1
                 
                 %f3 = ([File.pathname, '/Results','/', 'StressResults_finite_unconstrained', '.dat']);
                %dlmwrite(f3,finalResults,'precision',6);
                f3 = (['C:\code\Trac_MSM_Code\Data\', 'Results','\', 'StressResults_finite_unconstrained', '.dat']);
                dlmwrite(f3,finalResults,'precision',6);
                
             else
                 
                 %f3 = ([File.pathname, '/Results','/', 'StressResults_constrained', '.dat']);
                %dlmwrite(f3,finalResults,'precision',6);
                f3 = (['C:\code\Trac_MSM_Code\Data\', 'Results','\', 'StressResults_constrained', '.dat']);
                dlmwrite(f3,finalResults,'precision',6);
             end
             
            %save stress matrix 
            if (strcmpi(Settings.EvenStresses,'Yes')) == 1  
                  
                %save([File.pathname, '/Results','/', 'NewStresses', '.mat'], '-struct', 'NewStresses'); 
                save(['C:\code\Trac_MSM_Code\Data', '\Results','\', 'NewStresses', '.mat'], '-struct', 'NewStresses');
            end
             
             
         end
   
   
   
   close all
end
end