test1.m

clear; 
close all;
clc;

binout_filename='../LS-DYNA-sample/binout';
binin = get_binout_data(binout_filename);

%[figParents,figSelect]=struct2graph(binin,"All"); %graphically see contents of binin


%%
set(0,'defaulttextinterpreter','latex','defaulttextfontsize',7,'defaulttexthorizontalalignment','left','defaulttextverticalalignment','middle','defaulttextunits','data');
set(0,'defaultlegendinterpreter','latex');
%% matsum

Pids=binin.matsum.metadata.ids; %ids of parts
matsum=binin.matsum.data;  %rmfield(binin,'matsum');
t=matsum.time; % time vector for matsum database
ke=matsum.kinetic_energy; % kinetic energy time histories of part 1 (column 1) and part 2 (column 2)
ie=matsum.internal_energy; % internal energy
hge=matsum.hourglass_energy; % hourglass energy


figure(1); clf; cla; 
fig=gcf; ax=gca;
hold(ax,"on");

plot(ax,t,ke);
txt=compose('KE Part %d',Pids);
plot(ax,t,ie);
txt(end+1:end+length(Pids))=compose('IE Part %d',Pids);

plot(t,hge);
txt(end+1:end+length(Pids))=compose('HGE Part %d',Pids);

% plot(binin.glstat.data.time,binin.glstat.data.total_energy,'-k'); % total energy
% txt(end+1)={'TE'};

hold(ax,'off');
legend(ax,txt,'Location','northeast');
xlabel(ax,'Time [s]'); 
ylabel(ax,'Energy [N.m]');
set(ax.XAxis,'Exponent',-3,'TickLabelFormat','%.2f');
set(ax.YAxis,'TickLabelFormat','%.1f');
ax=tidyAxes(ax); % local fun, given at the bottom of this script

folderName='../figs';
if exist(folderName,'dir')==0
    mkdir(folderName);
end
figFileName=[folderName,'/','energies'];
fig=printFig(fig,figFileName,["pdf","svg"]); % another local fun, given at the bottom of this script



%% nodout

Nids=binin.nodout.metadata.ids; % ids of nodes
nodout=binin.nodout.data; %rmfield(binin,'nodout');
x0=nodout.x_coordinate(1,:)'; % inital x-coordinate
y0=nodout.y_coordinate(1,:)'; % ... y-coordinate

% search for node nearest to (xq,yq)
xq=0.5;
yq=0.2;
[~,idq]=min(sqrt((x0-xq).^2+(y0-yq).^2));
idq=idq(1); 
xq=x0(idq); yq=y0(idq); % actual coordinate of node
nid=Nids(idq); % actual id of node


t=nodout.time; % time vector of nodout database
uy=nodout.y_displacement(:,idq); % y-displacement time-history at nid
vy=nodout.y_velocity(:,idq); % y-velocity ...

figure(2); clf; cla;
fig=gcf; ax=gca;
plot(t,uy);
txt=compose('$u_y$ Node %d',nid);

legend(ax,txt);
xlabel(ax,'Time [s]');
ylabel(ax,'y-Displacement [m]');
set(ax.XAxis,'Exponent',-3,'TickLabelFormat','%.2f');
set(ax.YAxis,'Exponent',-3,'TickLabelFormat','%.2f');
ax=tidyAxes(ax);




%% displaying the model geometry and motion

connec=binin.control.connec.shell; %element-node connectivity array for shell elems 
x=nodout.x_coordinate; % 2D array of x-coordinate time histories of all nodes
y=nodout.y_coordinate; % ... y-coordinate ...
faces=connec(1:4,:)'; % for patch
vertices=[x(1,:)',y(1,:)']; % for patch 

vx=nodout.x_velocity;
vy=nodout.y_velocity;
vResultant=sqrt(vx.^2+vy.^2); % resultant velocity
vResultant(abs(vResultant)<=eps)=0;


figure(3); clf; cla;
fig=gcf; 
tl=tiledlayout(fig,1,1,'TileSpacing','tight','Padding','tight');
ax=nexttile;
axis(ax,"equal");
xlim([min(x,[],"all"),max(x,[],"all")]);
ylim([min(y,[],"all"),max(y,[],"all")]);
fig.Units='centimeters';
ax.Units='centimeters'; 
axPos=ax.tightPosition(IncludeLabels=true);
fig.Position(3:4)=axPos(3:4);



xlabel(ax,'x-coordinate [m]');
ylabel(ax,'y-coordinate [m]');
ax.TickDir="both";
ax.Color='none';
ax.set('FontSize',8);

pa =patch('Faces',faces,'Vertices',vertices,'FaceColor','interp','EdgeColor','none');
colorVar='Resultant Velocity';

pa.FaceVertexCData=vResultant(1,:)';
colormap(ax,'turbo');
cb=colorbar(ax); cb.Label.set('String',[colorVar,' [m/s]'],'Interpreter','none');
cb.Layout.Tile='east';
clim(ax,[min(vResultant,[],"all"),max(vResultant,[],"all")]);
txt=compose('Time = %7.6f [s]',t(1));
te=text(0.02,1,txt,'FontSize',8,'VerticalAlignment','top');

ax.NextPlot="replacechildren";
fig.ToolBar="none";
F=struct('cdata',[],'colormap',[]);
frameCount=min([120,length(t)]);
F(frameCount)=struct('cdata',[],'colormap',[]);
Im=cell(1,frameCount);

for i=1:frameCount
    pa.Vertices=[x(i,:)',y(i,:)'];
    pa.FaceVertexCData=vResultant(i,:)';
    te.String=compose('Time = %7.6f [s]',t(i));
    F(i)=getframe(fig);
    Im{i}=frame2im(F(i));
    pause(1/90);
end
F=[F([1 1 1]),F(2:end)];
frameAspectRatio=size(F(1).cdata,1)/size(F(1).cdata,2);


%% saving videos and gifs

figure(4); clf; cla;
fig=gcf; delete(gca());
fig.Position(4)=fig.Position(3)*(frameAspectRatio);
movie(fig,F,2,20,[0 0 0 0]);
videosFolder='../videos';
if exist(videosFolder,'dir')==0
    mkdir(videosFolder);
end
videoFileName=[videosFolder,'/','impact_color_resultant_vel'];
vid=VideoWriter([videoFileName,'.mp4'],'MPEG-4');
vid.FrameRate=12;
vid.Quality=90;
pads=2-mod(size(F(1).cdata),2);
pads(mod(size(F(1).cdata),2)==0)=0;
padW=pads(1); padH=pads(2);
if any([padW,padH]~=0)
    for i=1:length(F)
        F(i).cdata(end+1:end+padW,:,:)=0;
        F(i).cdata(:,end+1:end+padH,:)=0;
    end
end
open(vid);
writeVideo(vid,F);
close(vid);

for i=1:length(Im)
    [A,map]=rgb2ind(Im{i},256,'nodither');
    if i==1
        imwrite(A,map,[videoFileName,'.gif'],'LoopCount',Inf,'DelayTime',0.1);
    else
        imwrite(A,map,[videoFileName,'.gif'],'WriteMode',"append",'DelayTime',0.08);
    end
end



%% elout

elout=binin.elout.shell.data;

t2=elout.time; % time vector of elout database
if all(elout.nip(1,:)==1)
    % 1 integration point per elem. 
    pres=-(elout.sig_xx+elout.sig_yy+elout.sig_zz)/3; % pressure (+ve: compression)
else
    % integration points per element >1
    ipt=1; % chosen int. pt. 
    iptIdx=(1:size(elout.ids,2)).*elout.nip(1,:)-(elout.nip(1,:)-ipt); % indices of int. pt. #(ipt) for all elems
    pres=-(elout.sig_xx(:,iptIdx)+elout.sig_yy(:,iptIdx)+elout.sig_zz(:,iptIdx))/3; % pressure (+ve: compression)
    fprintf(1,'Note: there are %d int. pts. per elem. Chose int pt. #%d.\n',max(elout.nip(1,:)),ipt);
end
nid=4113; %node at mid point
[~,eids]=find(connec(1:4,:)==nid); % elem ids sharing node nid
presq=mean(pres(:,eids),2); % taking node pressure as average of those of surrounding elems

figure(5); clf; cla;
fig=gcf;  ax=gca;
yyaxis(ax,"left");
plot(ax,t2,presq)
txt=compose('P Node %d',nid);
text(ax,0.01e-3,0,{'$\uparrow$ compression'},'VerticalAlignment','bottom');
text(ax,0.01e-3,0,{'$\downarrow$ tension'},'VerticalAlignment','top');
ylabel(ax,'Pressure [Pa]');
set(ax.YAxis,'Exponent',9,'TickLabelFormat','%.1f');

yyaxis(ax,"right");
ydata=-nodout.y_velocity(:,nid);
plot(ax,t,ydata);
txt(end+1)=compose('$-\\dot{u}_y$ Node %d',nid);
text(ax,0.01e-3,0,{'$\uparrow$ down'},'VerticalAlignment','bottom');
text(ax,0.01e-3,0,{'$\downarrow$ up'},'VerticalAlignment','top');
ylabel(ax,'Velocity [m/s]');
set(ax.YAxis,'TickLabelFormat','%.1f');

legend(ax,txt);
xlabel(ax,'Time [s]');
set(ax.XAxis,'Exponent',-3,'TickLabelFormat','%.2f');
ax=tidyAxes(ax);

figFileName=[folderName,'/','pressure_and_vel_at_mid_node'];
fig=printFig(fig,figFileName,["pdf","svg"]);




%% displaying and annotating the input model

P3Nodes=connec(1:4,connec(end,:)==find(Pids==3));
P4Nodes=connec(1:4,connec(end,:)==find(Pids==4));
P3LB=[min(x0(P3Nodes(:))),min(y0(P3Nodes(:)))];
P4LB=[min(x0(P4Nodes(:))),min(y0(P4Nodes(:)))];


figure(6); clf; cla;
fig=gcf; ax=gca;
fig.Units='centimeters'; fig.Position(3:4)=[14,10];
axis(ax,"equal");
xlim([min(x,[],"all"),max(x,[],"all")]);
ylim([min(y,[],"all"),max(y,[],"all")]);

xlabel(ax,'x-coordinate [m]');
ylabel(ax,'y-coordinate [m]');
ax.TickDir="both";
ax.Color='none';
ax.set('FontSize',8);

faces=connec(1:4,:)'; % for patch
vertices=[x(1,:)',y(1,:)']; % for patch 

pa =patch('Faces',faces,'Vertices',vertices,'FaceColor','flat','EdgeColor','k');
pa.FaceVertexCData=10*connec(end,:)'; % part id indices
pa.FaceAlpha=0.8; pa.EdgeAlpha=0.5; pa.LineWidth=0.2;
colormap(ax,jet(4));

hold(ax,'on');
for mrk=["o","+"]
plot(ax,xq,yq,mrk,'MarkerFaceColor','none','MarkerEdgeColor','w','MarkerSize',6);
end
hold(ax,'off');     
txt=compose(' Node %d',nid);
fontsz=10;

th=text(ax,xq+0.02,yq,txt,'Color','w','FontWeight','normal','FontSize',fontsz);
% adding parts labels
txt=compose('Part %d',Pids);
text(ax,[P3LB(1);P4LB(1)]+0.01,[P3LB(2);P4LB(2)]+0.01,txt,'FontSize',fontsz,'Color','w','VerticalAlignment','bottom','Margin',4);
ax.Units='centimeters';
fig.Units='centimeters';
axPos=ax.tightPosition(IncludeLabels=true);
fig.Position(3:4)=axPos(3:4);
ax.Position(1:2)=ax.Position(1:2)+([0 0]-axPos(1:2));

txtID=findobj(ax.Children,'Type','Text');
for txtid=txtID'

tPos=txtid.Extent'+[0 0 txtid.Extent(1:2)]'+0.002*[-1 -0.25 1 0.25]';
tPos=reshape(tPos,2,[]);
[tX,tY]=ndgrid(tPos(1,:),tPos(2,:)); tPos=[tX(:),tY(:)];
txtpa=patch(ax,'Faces',[1,2,4,3],'Vertices',tPos,'FaceColor','k','FaceAlpha',0.3,'EdgeColor','k','EdgeAlpha',1,'LineWidth',0.2);
end
uistack(txtID,'top');

figFileName=[folderName,'/','initial_geometry_with_labels'];
fig=printFig(fig,figFileName,["pdf","svg"]);




%%
yq=[0,0.2,0.4];
xq=repmat(0.5,size(yq));

[~,idq]=min(sqrt((x0-xq).^2+(y0-yq).^2));
% idq=idq(1); 
xq=x0(idq); yq=y0(idq); % actual coordinate of node
nid=Nids(idq); % actual id of node
nid=double(nid);
pid=1;
[~,eidx]=ismember(connec(1:4,connec(end,:)==pid),nid); % elem ids sharing node nid
pres1=pres(:,connec(end,:)==pid);
presq=zeros(length(t),length(nid));
for i=1:length(nid)
    [~,eids]=find(eidx==i);
    presq(:,i)=mean(pres1(:,eids),2); % Note: for surface nodes, this is without extrapolation, i.e. pressures = nearest integration point vertically.
end

figure(7); clf; cla;
fig=gcf; ax=gca;
plot(ax,t,presq);
txt=compose('P(%g,%g), Node %d',[xq,yq,nid']);
legend(ax,txt,'Location','northeast');
xlabel(ax,'Time [s]'); ylabel(ax,'Pressure [Pa]');
set([ax.XAxis],'Exponent',-3,'TickLabelFormat','%.2f');
set(ax.YAxis,'Exponent',9,'TickLabelFormat','%.1f');
ax=tidyAxes(ax);
figFileName=[folderName,'/','pressure_near_edges_and_mid'];
fig=printFig(fig,figFileName,["pdf","svg"]);







%% local fun's

function [axOut]=tidyAxes(axIn)
% re-set defaults of axes:
set(axIn,'XGrid','on','YGrid','on','GridAlpha',0.12);
set(axIn,'XMinorGrid','on','YMinorGrid','on','MinorGridAlpha',0.11);
set(axIn,'XMinorTick','on','YMinorTick','on','TickDir','in');
set(axIn,'FontName','Times','FontSize',10);
set([axIn.XAxis.Label,axIn.YAxis.Label],'Interpreter','latex');
set(axIn,'Box','on');
set([allchild(axIn)],'LineWidth',0.75);
%
% remove un-needed white space (margins):
posTight=tightPosition(axIn,IncludeLabels=true);
axIn.Position(1:2)=axIn.Position(1:2)+([0 0]-posTight(1:2));
axIn.Position(3:4)=axIn.Position(3:4)+([1 1]-posTight(3:4));
%
fig=axIn.Parent;
leg=findobj(fig.Children,'Type','Legend');
if ~isempty(leg) && any(string(leg.Location)==["northeast","northwest","southeast","southwest"])
leg.Units='normalized'; leg.Interpreter='latex';
axIn.Units='normalized';
gap(1)=sum(axIn.Position([1 3]))-sum(leg.Position([1 3])); %right
gap(2)=sum(axIn.Position([2 4]))-sum(leg.Position([2 4])); %top
gap(3)=axIn.Position([1])-leg.Position([1]); %left
gap(4)=axIn.Position([2])-leg.Position([2]); %bottom
[~,gminId]=mink(abs(gap),2);
leg.Position(1:2)=leg.Position(1:2)+0.4*gap(gminId);
end
axOut=axIn;
end



function [fig]=printFig(fig,filename,fmt)
if nargin<3 || (nargin==3 && isempty(fmt))
    fmt="pdf";
end
fmt=string(fmt);
filename=erase(filename,regexpPattern('\.[a-zA-Z0-9]+$'));
filename=char(filename);
%
set(fig,'InvertHardCopy','off');
set(fig,'Color','w');
if any(fmt=="pdf")
    fig.Units="centimeters";
    fig.PaperUnits="centimeters";
    fig.PaperPositionMode="manual";
    fig.PaperSize=fig.Position(3:4); %make page size = fig size
    fig.PaperPosition=[0 0 fig.Position(3:4)]; % zero margins
    %will plot later below
end
%
% scalable vector graphics (best vectorised fig for web)  
if any(fmt=="svg")
    figFileName=filename;
    print(fig,figFileName,'-dsvg');
end
%
% raster graphics
for fmt0=fmt(ismember(fmt,["png","jpeg","jpg"]))
    figFileName=[filename,'.',char(fmt0)];
    exportgraphics(fig,figFileName,'ContentType','image','Resolution',600);
end
%
% vector graphics
for fmt0=fmt(ismember(fmt,["pdf","emf","eps"]))
    figFileName=[filename,'.',char(fmt0)];
    exportgraphics(fig,figFileName,'ContentType','vector');
end
set(fig,'Color','factory');
end