clc
close all
clear all
%*************************** NOTE ******************************
% 1) The main file for cell angle 15 degrees is MAIN_Motion_RepresentativeVolume_Theta15.m
% 2) The following function files are required to execute the main file
    % (a) draw.m                : Function to visualise the undeformed configuration of
    %                             the structure being simulated. The representative
    %                             system of cell walls are highlighted in black.
    
    % (b) animation.m           : Function to visualise the deformation process.
    
    % (c) find_AreaCell         : Function to compute the area of the cell
    %                             based on the coordinates of the vertices of the cell
    
    % (d) gas_constitutive_reln : Function to compute the gas pressure from
    %                             the constitutive relation of the gas which has been
    %                             described in the paper 
    
% 3) Before running the main file, the limits of the graph needs to checked
%    in the function files 'draw.m', and 'animation.m'. 
% 4) By default the function calls for 'draw.m' and 'animation.m' are
%    commented out. Uncomment the function calls on lines 142 and 232 for visualising the
%    results.
%****************************************************************
%%
%********** Enter the geometrical parameters of the honeycomb ************
l = 3.59; % length of mid-section of the inclined walls
h = 5.23; % length of mid-section of the vertical walls
thickness = 0.344; % cell wall thickness
theta = 15*pi/180; % cell angle

L_X1 = 76.632847; % overall length of the structure along X1 direction
%%
%****** M and c values extracted for the honeycomb with 11x11 cells ******

M = 1.056074; % M is the slope of the Vs-Vi relation
c = 14824; % c is the y-intercept in mm/s

%*** Calculate shock velocity based on the strain-rate/impact velocity ***
for i=1:1:13
    Strain_rate(i) = 400 + (i-1)*100;
    V_impact_array(i) = L_X1*Strain_rate(i);
    Shock_Velocity(i) = M*V_impact_array(i)+c; 
end
%%
index_SR = 5;
V_impact = V_impact_array(index_SR); % impact velocity corresponding to chosen index of strainrate array
Vs = Shock_Velocity(index_SR); % shock velocity corresponding to chosen index of strainrate array
%%
%*************** size of the structure being simulated ***************
rows = 10;
columns = 4;
num_nodes = 40;
%%
%************** initialising the x and y coordinates of the nodes *********
counter_nodes = 1;
for i=1:1:10
    for j=1:1:4
        Level(i,j) = counter_nodes;
        if i == 1
            Struct_NodalData(counter_nodes).x = 0 + (j-1)*2*l*cos(theta);
            Struct_NodalData(counter_nodes).y = 0;
        elseif i==2 || i==3
            Struct_NodalData(counter_nodes).x = l*cos(theta) + (j-1)*2*l*cos(theta);
            if i==2
                Struct_NodalData(counter_nodes).y = l*sin(theta);
            else
                Struct_NodalData(counter_nodes).y = l*sin(theta) + h;
            end
        elseif i==4 || i==5
            Struct_NodalData(counter_nodes).x = 2*l*cos(theta) + (j-1)*2*l*cos(theta);
            if i==4
                Struct_NodalData(counter_nodes).y = 2*l*sin(theta) + h;
            else
                Struct_NodalData(counter_nodes).y = 2*(l*sin(theta) + h);
            end
        elseif i==6 || i==7
            Struct_NodalData(counter_nodes).x = 3*l*cos(theta) + (j-1)*2*l*cos(theta);
            if i==6
                Struct_NodalData(counter_nodes).y = 3*l*sin(theta) + 2*h;
            else
                Struct_NodalData(counter_nodes).y = 3*(l*sin(theta) + h);
            end
        elseif i==8
            Struct_NodalData(counter_nodes).x = 4*l*cos(theta) + (j-1)*2*l*cos(theta);
            Struct_NodalData(counter_nodes).y = 4*l*sin(theta) + 3*h;
            
        elseif i==9
            Struct_NodalData(counter_nodes).x = 0 + (j-1)*2*l*cos(theta);
            Struct_NodalData(counter_nodes).y = -h;
        else
            Struct_NodalData(counter_nodes).x = l*cos(theta) + (j-1)*2*l*cos(theta);
            Struct_NodalData(counter_nodes).y = -h - l*sin(theta);
            
        end
        counter_nodes = counter_nodes + 1;
    end
end
%%
%************* Defining the connectivity between cell walls **************
counter_elements = 1;
for i=1:1:7
    for j=1:1:columns
        Connectivity(counter_elements).a = Level(i,j);
        Connectivity(counter_elements).b = Level(i+1,j);
        counter_elements = counter_elements + 1;
    end
end

for i=2:2:8
    for j=1:1:columns-1
        Connectivity(counter_elements).a = Level(i,j);
        Connectivity(counter_elements).b = Level(i-1,j+1);
        counter_elements = counter_elements + 1;
    end
    
end

i=9;
for j=1:1:columns
    Connectivity(counter_elements).a = Level(1,j);
    Connectivity(counter_elements).b = Level(i,j);
    counter_elements = counter_elements + 1;
end

i=10;
for j=1:1:columns
    Connectivity(counter_elements).a = Level(9,j);
    Connectivity(counter_elements).b = Level(i,j);
    counter_elements = counter_elements + 1;
end

for j=1:1:columns-1
    Connectivity(counter_elements).a = Level(i,j);
    Connectivity(counter_elements).b = Level(9,j+1);
    counter_elements = counter_elements + 1;
end
%%
% ******** function to visualise the undeformed configuration ************
%                    of the structure being simulated
draw(Struct_NodalData,Connectivity);
%%
%******************* Defining the time steps ******************
total_time = 10.2*l*cos(theta)/Vs;
delta_t = l*cos(theta)/(Vs*12);  % The time taken by the shock front  
                                % to travel the distance between succesive
                                % columns is discretised into 12 time
                                % intervals

num_intervals = round(total_time/delta_t);  % size of the time array
time(1) = 0; 
for i=2:1:num_intervals
    time(i) = time(i-1) + delta_t;
end
%%
%****** Calculating the evolution of X and Y coordinates of various ******
%          nodes based on the shock velocity and their positions 
for i=1:2:rows-1
    for j=1:1:columns
        node_number = Level(i,j);
        delay_shock(i,j) = (Struct_NodalData(node_number).x - 0)/Vs;
        index_delay_shock(i,j) = find(round(time,9) == round(delay_shock(i,j),9));
        
        Velocity_X(node_number,1:index_delay_shock(i,j)) = 0;
        
        for k=index_delay_shock(i,j):1:num_intervals
            
            Velocity_X(node_number,k) = V_impact;
        end
    end
end

for i=2:2:rows
    for j=1:1:columns
        node_number = Level(i,j);
        delay_shock(i,j) = (Struct_NodalData(node_number).x - 0)/Vs;
        index_delay_shock(i,j) = find(round(time,9) == round(delay_shock(i,j),9));
        
        previous_level_adjacent_node_number = Level(i-1,j);
        delay_shock_adjacentnode(i,j) = (Struct_NodalData(previous_level_adjacent_node_number).x - 0)/Vs;
        index_delay_shock_adjacentnode(i,j) = find(round(time,9) == round(delay_shock_adjacentnode(i,j),9));
        
        Velocity_X(node_number,1:index_delay_shock_adjacentnode(i,j)) = 0;
        Velocity_X(node_number,index_delay_shock_adjacentnode(i,j):index_delay_shock(i,j)) = V_impact/2;
        
        for k=index_delay_shock(i,j):1:num_intervals
            
            Velocity_X(node_number,k) = V_impact;
        end
    end
end
%%
for i=1:1:rows
    for j=1:1:columns
        counter = 1;
        node_number = Level(i,j);
        for k=1:1:num_intervals
            if k == 1
                X_Coordinate(node_number,k) = Struct_NodalData(node_number).x;
                Y_Coordinate(node_number,k) = Struct_NodalData(node_number).y;
            else
                X_Coordinate(node_number,k) = X_Coordinate(node_number,k-1) + Velocity_X(node_number,k-1)*delta_t;
                if mod(i,2) ~=0
                    Y_Coordinate(node_number,k) = Y_Coordinate(node_number,1);
                else
                    previous_level_adjacent_node_number = Level(i-1,j);
                    if k <= index_delay_shock_adjacentnode(i,j)
                        Y_Coordinate(node_number,k) = Y_Coordinate(node_number,1);
                    elseif k > index_delay_shock_adjacentnode(i,j) && k <= index_delay_shock(i,j)
                        %                         i
                        %                         j
                        %                         X_Coordinate(Level(i-1,j),1) + 2*l*cos(theta)
                        %                         X_max = 2*j*l*cos(theta)
                        X_max = X_Coordinate(Level(i-1,j),1) + 2*l*cos(theta);
                        if i ~= 10
                            Y_Coordinate(node_number,k) = Y_Coordinate(previous_level_adjacent_node_number,1) + sqrt(l^2-(X_max-X_Coordinate(previous_level_adjacent_node_number,k))^(2)*0.25);
                        else
                            Y_Coordinate(node_number,k) = Y_Coordinate(previous_level_adjacent_node_number,1) - sqrt(l^2-(X_max-X_Coordinate(previous_level_adjacent_node_number,k))^(2)*0.25);
                        end
                    else
                        Y_Coordinate(node_number,k) = Y_Coordinate(node_number,index_delay_shock(i,j));
                    end
                end
            end
        end
    end
end
%%
%*** function to visualise the simplified deformation process with time ***
clc
animation(X_Coordinate,Y_Coordinate,Connectivity,Struct_NodalData);
clc
%%
%*** Each cell has 6 vertices with the respective node numbers assigned ***
Cell_Number(1,1) = 9;
Cell_Number(1,2) = 13;
Cell_Number(1,3) = 10;
Cell_Number(1,4) = 6;
Cell_Number(1,5) = 2;
Cell_Number(1,6) = 5;

Cell_Number(2,1) = 10;
Cell_Number(2,2) = 14;
Cell_Number(2,3) = 11;
Cell_Number(2,4) = 7;
Cell_Number(2,5) = 3;
Cell_Number(2,6) = 6;

Cell_Number(3,1) = 11;
Cell_Number(3,2) = 15;
Cell_Number(3,3) = 12;
Cell_Number(3,4) = 8;
Cell_Number(3,5) = 4;
Cell_Number(3,6) = 7;

Cell_Number(4,1) = 17;
Cell_Number(4,2) = 21;
Cell_Number(4,3) = 18;
Cell_Number(4,4) = 14;
Cell_Number(4,5) = 10;
Cell_Number(4,6) = 13;

Cell_Number(5,1) = 18;
Cell_Number(5,2) = 22;
Cell_Number(5,3) = 19;
Cell_Number(5,4) = 15;
Cell_Number(5,5) = 11;
Cell_Number(5,6) = 14;

Cell_Number(6,1) = 19;
Cell_Number(6,2) = 23;
Cell_Number(6,3) = 20;
Cell_Number(6,4) = 16;
Cell_Number(6,5) = 12;
Cell_Number(6,6) = 15;

Cell_Number(7,1) = 25;
Cell_Number(7,2) = 29;
Cell_Number(7,3) = 26;
Cell_Number(7,4) = 22;
Cell_Number(7,5) = 18;
Cell_Number(7,6) = 21;

Cell_Number(8,1) = 26;
Cell_Number(8,2) = 30;
Cell_Number(8,3) = 27;
Cell_Number(8,4) = 23;
Cell_Number(8,5) = 19;
Cell_Number(8,6) = 22;

Cell_Number(9,1) = 27;
Cell_Number(9,2) = 31;
Cell_Number(9,3) = 28;
Cell_Number(9,4) = 24;
Cell_Number(9,5) = 20;
Cell_Number(9,6) = 23;

Cell_Number(10,1) = 1;
Cell_Number(10,2) = 5;
Cell_Number(10,3) = 2;
Cell_Number(10,4) = 34;
Cell_Number(10,5) = 37;
Cell_Number(10,6) = 33;

Cell_Number(11,1) = 2;
Cell_Number(11,2) = 6;
Cell_Number(11,3) = 3;
Cell_Number(11,4) = 35;
Cell_Number(11,5) = 38;
Cell_Number(11,6) = 34;

Cell_Number(12,1) = 3;
Cell_Number(12,2) = 7;
Cell_Number(12,3) = 4;
Cell_Number(12,4) = 36;
Cell_Number(12,5) = 39;
Cell_Number(12,6) = 35;
%%
% The area between the mid-sections and inner sections need to subtracted
%         to calculate the vlume occupied by the entrappad gas
%*********** lengths of the inner sections ***********
l_inner = l-(thickness/2)*tan(15*pi/180)-(thickness/2)*(1/tan(52.5*pi/180));
h_inner = h-(thickness/2)*(2/tan(52.5*pi/180));

subtract_area = 0.5*(4*l+2*h+4*l_inner+2*h_inner)*0.5*thickness;
Volume_Cells = find_AreaCell(Cell_Number,X_Coordinate,Y_Coordinate,subtract_area);
num_cells = size(Volume_Cells,1);
%%
%******** Computing the gas pressure from the constitutive relation *******
%                           of the ideal gas 
clc
counter_10_11_12  = 1;
for i=1:1:num_cells
    for j=1:1:num_intervals-1
        if round(Volume_Cells(i,j),12) ~= round(Volume_Cells(i,j+1),12)
            StartEnd_Index_PressureCalculation(i,1) = j;
            break;
        end
    end
    
    for k=num_intervals:-1:2
        if round(Volume_Cells(i,k),12) ~= round(Volume_Cells(i,k-1),12)
            StartEnd_Index_PressureCalculation(i,2) = k;
            break;
        end
    end
    
    if i<=9
        V(i,:) = Volume_Cells(i,StartEnd_Index_PressureCalculation(i,1):StartEnd_Index_PressureCalculation(i,2));
    else
        V_10_11_12(counter_10_11_12,:) = Volume_Cells(i,StartEnd_Index_PressureCalculation(i,1):StartEnd_Index_PressureCalculation(i,2));
        counter_10_11_12 = counter_10_11_12 + 1;
    end
end

clc
for i=1:1:9
    V0 = Volume_Cells(i,StartEnd_Index_PressureCalculation(i,1));
    P0 = 0.2*10^6;
    
    if index_SR<=7
        
        for k=1:1:size(V,2)
            if V(i,k) < V(i,k+1)
                Break_point(i) = k;
                break;
            end
        end
        
        V_part1(i,:) = V(i,1:Break_point(i));
        V_part2(i,:) = cat(2,V_part1(i,1),V(i,Break_point(i)+1:end));
        
        V1 = V_part1(i,:);
        V2 = V_part2(i,:);
        
        [V1,P1] = ode23(@(V1,P1)gas_constitutive_reln(V1,P1),V1,P0);
        for j=1:1:size(P1,1)
            Pressure_Cells_123456(i,j) = P1(j);
        end
        
        [V2,P2] = ode23(@(V2,P2)gas_constitutive_reln(V2,P2),V2,P0);
        temp = P2(2:end);
        for j=1:1:size(temp,1)
            Pressure_Cells_123456(i,j+size(P1,1))= temp(j);
        end
        
        Pressure_Cells(i,1:StartEnd_Index_PressureCalculation(i,1)-1) = Pressure_Cells_123456(i,1);
        Pressure_Cells(i,StartEnd_Index_PressureCalculation(i,1):StartEnd_Index_PressureCalculation(i,2)) = Pressure_Cells_123456(i,:);
        Pressure_Cells(i,StartEnd_Index_PressureCalculation(i,2):num_intervals) = Pressure_Cells_123456(i,1);
    else
        V1 = V(i,:);
        [V1,P1] = ode23(@(V1,P1)gas_constitutive_reln(V1,P1),V1,P0);
        for j=1:1:size(P1,1)
            Pressure_Cells_1to8(i,j) = P1(j);
        end
        
        Pressure_Cells(i,1:StartEnd_Index_PressureCalculation(i,1)-1) = Pressure_Cells_1to8(i,1);
        Pressure_Cells(i,StartEnd_Index_PressureCalculation(i,1):StartEnd_Index_PressureCalculation(i,2)) = Pressure_Cells_1to8(i,:);
        Pressure_Cells(i,StartEnd_Index_PressureCalculation(i,2):num_intervals) = Pressure_Cells_1to8(i,end);
    end
end

P_History(1,:) = Pressure_Cells(4,:);
P_History(2,:) = Pressure_Cells(2,:);
P_History(3,:) = Pressure_Cells(7,:);
P_History(4,:) = Pressure_Cells(3,:);
P_History(5,:) = Pressure_Cells(8,:);
P_History(6,:) = Pressure_Cells(6,:);
P_History(7,:) = Pressure_Cells(5,:);
%%
% Time duration of crushing (Eqn 15 in paper) in the impulse-momentum theorem
Index_Initial_Time = index_delay_shock(7,1);
Index_Final_Time = index_delay_shock(4,2);
collapse_duration = time(Index_Final_Time) - time(Index_Initial_Time);

% Calculating the gauge pressure of cells for the above time duration

Cell_Pressure(1,:) = Pressure_Cells(4,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
Cell_Pressure(2,:) = Pressure_Cells(2,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
Cell_Pressure(3,:) = Pressure_Cells(7,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
Cell_Pressure(4,:) = Pressure_Cells(3,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
Cell_Pressure(5,:) = Pressure_Cells(8,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
Cell_Pressure(6,:) = Pressure_Cells(6,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
Cell_Pressure(7,:) = Pressure_Cells(5,Index_Initial_Time:Index_Final_Time)*10^-6 - 0.1;
%%
%*********** cell walls belonging to the representative system ***********
Required_elements = [10;14;35;21;33;18;11;36;22;15];
counter = 1;

for j=Index_Initial_Time:1:Index_Final_Time
    
    for i=1:1:size(Required_elements,1)
        
        node1 = Connectivity(Required_elements(i)).a;
        node2 = Connectivity(Required_elements(i)).b;
        
        Representative_CellWalls(counter).cellwall(i).XComp = X_Coordinate(node1,j) - X_Coordinate(node2,j);
        Representative_CellWalls(counter).cellwall(i).YComp = Y_Coordinate(node1,j) - Y_Coordinate(node2,j);
        Representative_CellWalls(counter).cellwall(i).UneditedLength = sqrt(Representative_CellWalls(counter).cellwall(i).XComp^2 + Representative_CellWalls(counter).cellwall(i).YComp^2);
        if (i==14 | (i==21) | (i==15) | (i==22))
            Representative_CellWalls(counter).cellwall(i).Length = Representative_CellWalls(counter).cellwall(i).UneditedLength - (thickness/tan(52.5*pi/180));
        else
            Representative_CellWalls(counter).cellwall(i).Length = Representative_CellWalls(counter).cellwall(i).UneditedLength - 0.5*thickness*(tan(15*pi/180)+(1/tan(52.5*pi/180)));
        end
        Representative_CellWalls(counter).cellwall(i).Angle = acos(Representative_CellWalls(counter).cellwall(i).YComp/Representative_CellWalls(counter).cellwall(i).UneditedLength)*180/pi;
        
    end
    counter = counter + 1;
end
%%
%**Computing the pressure force (see Eqn A.1 in the appendix of the paper**
num_frames_duringcollapse = size(Representative_CellWalls,2);
Required_Time = time(Index_Initial_Time:Index_Final_Time);

for i=1:1:num_frames_duringcollapse
    AB_integrand(i) = Cell_Pressure(2,i)*abs(cos(Representative_CellWalls(i).cellwall(1).Angle*pi/180)*Representative_CellWalls(i).cellwall(1).Length);
    BC_integrand(i) = Cell_Pressure(7,i)*abs(cos(Representative_CellWalls(i).cellwall(2).Angle*pi/180)*Representative_CellWalls(i).cellwall(2).Length);
    CD_integrand(i) = Cell_Pressure(3,i)*abs(cos(Representative_CellWalls(i).cellwall(3).Angle*pi/180)*Representative_CellWalls(i).cellwall(3).Length);
    DE_integrand(i) = Cell_Pressure(3,i)*abs(cos(Representative_CellWalls(i).cellwall(4).Angle*pi/180)*Representative_CellWalls(i).cellwall(4).Length);
    FB_integrand(i) = (Cell_Pressure(7,i)-Cell_Pressure(2,i))*abs(cos(Representative_CellWalls(i).cellwall(5).Angle*pi/180)*Representative_CellWalls(i).cellwall(5).Length);
    CG_integrand(i) = (Cell_Pressure(7,i)-Cell_Pressure(3,i))*abs(cos(Representative_CellWalls(i).cellwall(6).Angle*pi/180)*Representative_CellWalls(i).cellwall(6).Length);
    FI_integrand(i) = (Cell_Pressure(4,i)-Cell_Pressure(7,i))*abs(cos(Representative_CellWalls(i).cellwall(8).Angle*pi/180)*Representative_CellWalls(i).cellwall(8).Length);
    GJ_integrand(i) = (Cell_Pressure(5,i)-Cell_Pressure(7,i))*abs(cos(Representative_CellWalls(i).cellwall(10).Angle*pi/180)*Representative_CellWalls(i).cellwall(10).Length);
    GH_integrand(i) = (Cell_Pressure(5,i)-Cell_Pressure(3,i))*abs(cos(Representative_CellWalls(i).cellwall(7).Angle*pi/180)*Representative_CellWalls(i).cellwall(7).Length);
    IJ_integrand(i) = (Cell_Pressure(6,i)-Cell_Pressure(7,i))*abs(cos(Representative_CellWalls(i).cellwall(9).Angle*pi/180)*Representative_CellWalls(i).cellwall(9).Length);
end

Pressure_Force(1) = trapz(Required_Time,AB_integrand);
Pressure_Force(2) = trapz(Required_Time,BC_integrand);
Pressure_Force(3) = trapz(Required_Time,CD_integrand);
Pressure_Force(4) = trapz(Required_Time,DE_integrand);
Pressure_Force(5) = trapz(Required_Time,FB_integrand);
Pressure_Force(6) = trapz(Required_Time,CG_integrand);
Pressure_Force(7) = trapz(Required_Time,FI_integrand);
Pressure_Force(8) = trapz(Required_Time,GJ_integrand);
Pressure_Force(9) = trapz(Required_Time,GH_integrand);
Pressure_Force(10) = trapz(Required_Time,IJ_integrand);

PF = sum(Pressure_Force);
%%
%********** Calculating the stress given by Eqn. 24 in the paper **********
V = V_impact/1000;
delta_t_gas = collapse_duration;

term1_quasi = 30*((thickness/l)^2)/(2*sin(15*pi/180)*((h/l)+sin(15*pi/180)));
term2_inertial = (2*thickness*l + thickness*h)*2700*V*10^-9/(delta_t_gas*(2*l+h));
term3_PF = PF/(delta_t_gas*(2*l+h));

analytical_stress = term1_quasi + term2_inertial + term3_PF;

clc
fprintf('Cell angle = %.0f \n',theta*180/pi)
fprintf('Strain-rate = %.0f \n',Strain_rate(index_SR))
fprintf('Impact velocity = %f  mm/s\n',V_impact)
fprintf('Shock velocity = %f mm/s\n',Vs)
fprintf('Theoretical_stress = %f MPa\n',round(analytical_stress,4))
fprintf('Change index_SR on line 47 for different strain-rates\n')
