function displacement_processing close all load('compression_data.mat') % Read the compression data file % load('Cell_data.mat') % Read the cell data file % % Each data file contains: % nx, ny, nz, nt = no. grid points in x,y,z, directions and # time steps % xyz(1:3,nn) are xyz coords of grid point number nn. % uvw(1:3,nn,i) are the xyz components of displacemnt for point nn % at time step i % The order of the grid points is such that % xyz(i,1:nx) is the first row of points along the x direction, % with y=z=0 % xyz(i,nx+1:2*nx) is the second row of the top plane with z=0 % xyz(i,nx*ny+1) is the start of the second plane with z=dz % % Create figure figure1 = figure; for t = 1:nt % Loop over time steps for k = 1:nz for j = 1:ny for i = 1:nx nn = i+(j-1)*nx + (k-1)*ny*nx; x(i,j,k) = xyz(1,nn); y(i,j,k) = xyz(2,nn); z(i,j,k) = xyz(3,nn); u(i,j,k) = uvw(1,nn,t); % This plots x component of displ. end end end % % MATLAB slice expects the order of coords to be indexed % with y varying most rapidly, then x, then z. % P = [2 1 3]; X = permute(x, P); Y = permute(y, P); Z = permute(z, P); V = permute(u, P); clf % Clear the figure axes1 = axes('Parent',figure1,... 'Position',[0.13 0.14047619047619 0.648571428571429 0.78452380952381],... 'FontSize',14,... 'CLim',[-2.64384412765503 1.82826542854309]); view(axes1,[-37.5 30]); grid(axes1,'on'); hold(axes1,'all'); slice(X,Y,Z,V,xyz(2,end)/2,xyz(1,end)/2,xyz(3,end)/2) % Create xlabel xlabel('x (microns)','FontSize',14); % Create ylabel ylabel({'y (microns)'},'FontSize',14); % Create zlabel zlabel({'z (microns)'},'FontSize',14); % Create colorbar colorbar('peer',axes1,... [0.864177927927928 0.127272727272727 0.0225225225225226 0.690078870900789],... 'FontSize',14); % Create textbox annotation(figure1,'textbox',... [0.783494208494208 0.89821722113503 0.195945945945948 0.0639269406392694],... 'String',{'Displacements (microns)'},... 'HorizontalAlignment','center',... 'FontSize',14,... 'FitBoxToText','off',... 'LineStyle','none'); pause(0.1) end end % %================= SHAPE FUNCTIONS ================================== % % Calculates shape functions for various element types % function N = shapefunctions(nelnodes,ncoord,xi) ! nelnodes - no. nodes on the element ! ncoord - no. coords (2 for 2D, 3 for 3D) ! xi - coords of point where dNdxi is to be evaluated. N = zeros(nelnodes,1); % % 1D elements % if (ncoord == 1) if (nelnodes==2) N(1) = 0.5*(1.+xi(1)); N(2) = 0.5*(1.-xi(1)); elseif (nelnodes == 3) N(1) = -0.5*xi(1)*(1.-xi(1)); N(2) = 0.5*xi(1)*(1.+xi(1)); N(3) = (1.-xi(1))*(1.+xi(1)); end % % 2D elements % elseif (ncoord == 2) % % Triangular element % if ( nelnodes == 3 ) N(1) = xi(1); N(2) = xi(2); N(3) = 1.-xi(1)-xi(2); elseif ( nelnodes == 6 ) xi3 = 1.-xi(1)-xi(2); N(1) = (2.*xi(1)-1.)*xi(1); N(2) = (2.*xi(2)-1.)*xi(2); N(3) = (2.*xi3-1.)*xi3; N(4) = 4.*xi(1)*xi(2); N(5) = 4.*xi(2)*xi3; N(6) = 4.*xi3*xi(1); % % Rectangular element % elseif ( nelnodes == 4 ) N(1) = 0.25*(1.-xi(1))*(1.-xi(2)); N(2) = 0.25*(1.+xi(1))*(1.-xi(2)); N(3) = 0.25*(1.+xi(1))*(1.+xi(2)); N(4) = 0.25*(1.-xi(1))*(1.+xi(2)); elseif (nelnodes == 8) N(1) = -0.25*(1.-xi(1))*(1.-xi(2))*(1.+xi(1)+xi(2)); N(2) = 0.25*(1.+xi(1))*(1.-xi(2))*(xi(1)-xi(2)-1.); N(3) = 0.25*(1.+xi(1))*(1.+xi(2))*(xi(1)+xi(2)-1.); N(4) = 0.25*(1.-xi(1))*(1.+xi(2))*(xi(2)-xi(1)-1.); N(5) = 0.5*(1.-xi(1)*xi(1))*(1.-xi(2)); N(6) = 0.5*(1.+xi(1))*(1.-xi(2)*xi(2)); N(7) = 0.5*(1.-xi(1)*xi(1))*(1.+xi(2)); N(8) = 0.5*(1.-xi(1))*(1.-xi(2)*xi(2)); end % elseif (ncoord==3) if (nelnodes == 4) N(1) = xi(1); N(2) = xi(2); N(3) = xi(3); N(4) = 1.-xi(1)-xi(2)-xi(3); elseif (nelnodes == 10) xi4 = 1.-xi(1)-xi(2)-xi(3); N(1) = (2.*xi(1)-1.)*xi(1); N(2) = (2.*xi(2)-1.)*xi(2); N(3) = (2.*xi(3)-1.)*xi(3); N(4) = (2.*xi4-1.)*xi4; N(5) = 4.*xi(1)*xi(2); N(6) = 4.*xi(2)*xi(3); N(7) = 4.*xi(3)*xi(1); N(8) = 4.*xi(1)*xi4; N(9) = 4.*xi(2)*xi4; N(10) = 4.*xi(3)*xi4; elseif (nelnodes == 8) N(1) = (1.-xi(1))*(1.-xi(2))*(1.-xi(3))/8.; N(2) = (1.+xi(1))*(1.-xi(2))*(1.-xi(3))/8.; N(3) = (1.+xi(1))*(1.+xi(2))*(1.-xi(3))/8.; N(4) = (1.-xi(1))*(1.+xi(2))*(1.-xi(3))/8.; N(5) = (1.-xi(1))*(1.-xi(2))*(1.+xi(3))/8.; N(6) = (1.+xi(1))*(1.-xi(2))*(1.+xi(3))/8.; N(7) = (1.+xi(1))*(1.+xi(2))*(1.+xi(3))/8.; N(8) = (1.-xi(1))*(1.+xi(2))*(1.+xi(3))/8.; elseif (nelnodes == 20) N(1) = (1.-xi(1))*(1.-xi(2))*(1.-xi(3))*(-xi(1)-xi(2)-xi(3)-2.)/8.; N(2) = (1.+xi(1))*(1.-xi(2))*(1.-xi(3))*(xi(1)-xi(2)-xi(3)-2.)/8.; N(3) = (1.+xi(1))*(1.+xi(2))*(1.-xi(3))*(xi(1)+xi(2)-xi(3)-2.)/8.; N(4) = (1.-xi(1))*(1.+xi(2))*(1.-xi(3))*(-xi(1)+xi(2)-xi(3)-2.)/8.; N(5) = (1.-xi(1))*(1.-xi(2))*(1.+xi(3))*(-xi(1)-xi(2)+xi(3)-2.)/8.; N(6) = (1.+xi(1))*(1.-xi(2))*(1.+xi(3))*(xi(1)-xi(2)+xi(3)-2.)/8.; N(7) = (1.+xi(1))*(1.+xi(2))*(1.+xi(3))*(xi(1)+xi(2)+xi(3)-2.)/8.; N(8) = (1.-xi(1))*(1.+xi(2))*(1.+xi(3))*(-xi(1)+xi(2)+xi(3)-2.)/8.; N(9) = (1.-xi(1)^2)*(1.-xi(2))*(1.-xi(3))/4.; N(10) = (1.+xi(1))*(1.-xi(2)^2)*(1.-xi(3))/4.; N(11) = (1.-xi(1)^2)*(1.+xi(2))*(1.-xi(3))/4.; N(12) = (1.-xi(1))*(1.-xi(2)^2)*(1.-xi(3))/4.; N(13) = (1.-xi(1)^2)*(1.-xi(2))*(1.+xi(3))/4.; N(14) = (1.+xi(1))*(1.-xi(2)^2)*(1.+xi(3))/4.; N(15) = (1.-xi(1)^2)*(1.+xi(2))*(1.+xi(3))/4.; N(16) = (1.-xi(1))*(1.-xi(2)^2)*(1.+xi(3))/4.; N(17) = (1.-xi(1))*(1.-xi(2))*(1.-xi(3)^2)/4.; N(18) = (1.+xi(1))*(1.-xi(2))*(1.-xi(3)^2)/4.; N(19) = (1.+xi(1))*(1.+xi(2))*(1.-xi(3)^2)/4.; N(20) = (1.-xi(1))*(1.+xi(2))*(1.-xi(3)^2)/4.; end end end % %================= SHAPE FUNCTION DERIVATIVES ====================== % function dNdxi = shapefunctionderivs(nelnodes,ncoord,xi) ! nelnodes - no. nodes on the element ! ncoord - no. coords (2 for 2D, 3 for 3D) ! xi - coords of point where dNdxi is to be evaluated. dNdxi = zeros(nelnodes,ncoord); % % 1D elements % if (ncoord == 1) if (nelnodes==2) dNdxi(1,1) = 0.5; dNdxi(2,1) = -0.5; elseif (nelnodes == 3) dNdxi(1,1) = -0.5+xi(1); dNdxi(2,1) = 0.5+xi(1); dNdxi(3,1) = -2.*xi(1); end % % 2D elements % elseif (ncoord == 2) % % Triangular element % if ( nelnodes == 3 ) dNdxi(1,1) = 1.; dNdxi(2,2) = 1.; dNdxi(3,1) = -1.; dNdxi(3,2) = -1.; elseif ( nelnodes == 6 ) xi3 = 1.-xi(1)-xi(2); dNdxi(1,1) = 4.*xi(1)-1.; dNdxi(2,2) = 4.*xi(2)-1.; dNdxi(3,1) = -(4.*xi3-1.); dNdxi(3,2) = -(4.*xi3-1.); dNdxi(4,1) = 4.*xi(2); dNdxi(4,2) = 4.*xi(1); dNdxi(5,1) = -4.*xi(2); dNdxi(5,2) = -4.*xi(1); dNdxi(6,1) = 4.*xi3 - 4.*xi(1); dNdxi(6,2) = 4.*xi3 - 4.*xi(2); % % Rectangular element % elseif ( nelnodes == 4 ) dNdxi(1,1) = -0.25*(1.-xi(2)); dNdxi(1,2) = -0.25*(1.-xi(1)); dNdxi(2,1) = 0.25*(1.-xi(2)); dNdxi(2,2) = -0.25*(1.+xi(1)); dNdxi(3,1) = 0.25*(1.+xi(2)); dNdxi(3,2) = 0.25*(1.+xi(1)); dNdxi(4,1) = -0.25*(1.+xi(2)); dNdxi(4,2) = 0.25*(1.-xi(1)); elseif (nelnodes == 8) dNdxi(1,1) = 0.25*(1.-xi(2))*(2.*xi(1)+xi(2)); dNdxi(1,2) = 0.25*(1.-xi(1))*(xi(1)+2.*xi(2)); dNdxi(2,1) = 0.25*(1.-xi(2))*(2.*xi(1)-xi(2)); dNdxi(2,2) = 0.25*(1.+xi(1))*(2.*xi(2)-xi(1)); dNdxi(3,1) = 0.25*(1.+xi(2))*(2.*xi(1)+xi(2)); dNdxi(3,2) = 0.25*(1.+xi(1))*(2.*xi(2)+xi(1)); dNdxi(4,1) = 0.25*(1.+xi(2))*(2.*xi(1)-xi(2)); dNdxi(4,2) = 0.25*(1.-xi(1))*(2.*xi(2)-xi(1)); dNdxi(5,1) = -xi(1)*(1.-xi(2)); dNdxi(5,2) = -0.5*(1.-xi(1)*xi(1)); dNdxi(6,1) = 0.5*(1.-xi(2)*xi(2)); dNdxi(6,2) = -(1.+xi(1))*xi(2); dNdxi(7,1) = -xi(1)*(1.+xi(2)); dNdxi(7,2) = 0.5*(1.-xi(1)*xi(1)); dNdxi(8,1) = -0.5*(1.-xi(2)*xi(2)); dNdxi(8,2) = -(1.-xi(1))*xi(2); end % % 3D elements % elseif (ncoord==3) if (nelnodes == 4) dNdxi(1,1) = 1.; dNdxi(2,2) = 1.; dNdxi(3,3) = 1.; dNdxi(4,1) = -1.; dNdxi(4,2) = -1.; dNdxi(4,3) = -1.; elseif (nelnodes == 10) xi4 = 1.-xi(1)-xi(2)-xi(3); dNdxi(1,1) = (4.*xi(1)-1.); dNdxi(2,2) = (4.*xi(2)-1.); dNdxi(3,3) = (4.*xi(3)-1.); dNdxi(4,1) = -(4.*xi4-1.); dNdxi(4,2) = -(4.*xi4-1.); dNdxi(4,3) = -(4.*xi4-1.); dNdxi(5,1) = 4.*xi(2); dNdxi(5,2) = 4.*xi(1); dNdxi(6,2) = 4.*xi(3); dNdxi(6,3) = 4.*xi(2); dNdxi(7,1) = 4.*xi(3); dNdxi(7,3) = 4.*xi(1); dNdxi(8,1) = 4.*(xi4-xi(1)); dNdxi(8,2) = -4.*xi(1); dNdxi(8,3) = -4.*xi(1); dNdxi(9,1) = -4.*xi(2); dNdxi(9,2) = 4.*(xi4-xi(2)); dNdxi(9,3) = -4.*xi(2); dNdxi(10,1) = -4.*xi(3)*xi4; dNdxi(10,2) = -4.*xi(3); dNdxi(10,3) = 4.*(xi4-xi(3)); elseif (nelnodes == 8) dNdxi(1,1) = -(1.-xi(2))*(1.-xi(3))/8.; dNdxi(1,2) = -(1.-xi(1))*(1.-xi(3))/8.; dNdxi(1,3) = -(1.-xi(1))*(1.-xi(2))/8.; dNdxi(2,1) = (1.-xi(2))*(1.-xi(3))/8.; dNdxi(2,2) = -(1.+xi(1))*(1.-xi(3))/8.; dNdxi(2,3) = -(1.+xi(1))*(1.-xi(2))/8.; dNdxi(3,1) = (1.+xi(2))*(1.-xi(3))/8.; dNdxi(3,2) = (1.+xi(1))*(1.-xi(3))/8.; dNdxi(3,3) = -(1.+xi(1))*(1.+xi(2))/8.; dNdxi(4,1) = -(1.+xi(2))*(1.-xi(3))/8.; dNdxi(4,2) = (1.-xi(1))*(1.-xi(3))/8.; dNdxi(4,3) = -(1.-xi(1))*(1.+xi(2))/8.; dNdxi(5,1) = -(1.-xi(2))*(1.+xi(3))/8.; dNdxi(5,2) = -(1.-xi(1))*(1.+xi(3))/8.; dNdxi(5,3) = (1.-xi(1))*(1.-xi(2))/8.; dNdxi(6,1) = (1.-xi(2))*(1.+xi(3))/8.; dNdxi(6,2) = -(1.+xi(1))*(1.+xi(3))/8.; dNdxi(6,3) = (1.+xi(1))*(1.-xi(2))/8.; dNdxi(7,1) = (1.+xi(2))*(1.+xi(3))/8.; dNdxi(7,2) = (1.+xi(1))*(1.+xi(3))/8.; dNdxi(7,3) = (1.+xi(1))*(1.+xi(2))/8.; dNdxi(8,1) = -(1.+xi(2))*(1.+xi(3))/8.; dNdxi(8,2) = (1.-xi(1))*(1.+xi(3))/8.; dNdxi(8,3) = (1.-xi(1))*(1.+xi(2))/8.; elseif (nelnodes == 20) dNdxi(1,1) = (-(1.-xi(2))*(1.-xi(3))*(-xi(1)-xi(2)-xi(3)-2.)-(1.-xi(1))*(1.-xi(2))*(1.-xi(3)))/8.; dNdxi(1,2) = (-(1.-xi(1))*(1.-xi(3))*(-xi(1)-xi(2)-xi(3)-2.)-(1.-xi(1))*(1.-xi(2))*(1.-xi(3)))/8.; dNdxi(1,3) = (-(1.-xi(1))*(1.-xi(2))*(-xi(1)-xi(2)-xi(3)-2.)-(1.-xi(1))*(1.-xi(2))*(1.-xi(3)))/8.; dNdxi(2,1) = ((1.-xi(2))*(1.-xi(3))*(xi(1)-xi(2)-xi(3)-2.)+(1.+xi(1))*(1.-xi(2))*(1.-xi(3)))/8.; dNdxi(2,2) = (-(1.+xi(1))*(1.-xi(3))*(xi(1)-xi(2)-xi(3)-2.)-(1.+xi(1))*(1.-xi(2))*(1.-xi(3)))/8.; dNdxi(2,3) = (-(1.+xi(1))*(1.-xi(2))*(xi(1)-xi(2)-xi(3)-2.)-(1.+xi(1))*(1.-xi(2))*(1.-xi(3)))/8.; dNdxi(3,1) = ((1.+xi(2))*(1.-xi(3))*(xi(1)+xi(2)-xi(3)-2.)+(1.+xi(1))*(1.+xi(2))*(1.-xi(3)))/8.; dNdxi(3,2) = ((1.+xi(1))*(1.-xi(3))*(xi(1)+xi(2)-xi(3)-2.)+(1.+xi(1))*(1.+xi(2))*(1.-xi(3)))/8.; dNdxi(3,3) = (-(1.+xi(1))*(1.+xi(2))*(xi(1)+xi(2)-xi(3)-2.)-(1.+xi(1))*(1.+xi(2))*(1.-xi(3)))/8.; dNdxi(4,1) = (-(1.+xi(2))*(1.-xi(3))*(-xi(1)+xi(2)-xi(3)-2.)-(1.-xi(1))*(1.+xi(2))*(1.-xi(3)))/8.; dNdxi(4,2) = ((1.-xi(1))*(1.-xi(3))*(-xi(1)+xi(2)-xi(3)-2.)+(1.-xi(1))*(1.+xi(2))*(1.-xi(3)))/8.; dNdxi(4,3) = (-(1.-xi(1))*(1.+xi(2))*(-xi(1)+xi(2)-xi(3)-2.)-(1.-xi(1))*(1.+xi(2))*(1.-xi(3)))/8.; dNdxi(5,1) = (-(1.-xi(2))*(1.+xi(3))*(-xi(1)-xi(2)+xi(3)-2.)-(1.-xi(1))*(1.-xi(2))*(1.+xi(3)))/8.; dNdxi(5,2) = (-(1.-xi(1))*(1.+xi(3))*(-xi(1)-xi(2)+xi(3)-2.)-(1.-xi(1))*(1.-xi(2))*(1.+xi(3)))/8.; dNdxi(5,3) = ((1.-xi(1))*(1.-xi(2))*(-xi(1)-xi(2)+xi(3)-2.)+(1.-xi(1))*(1.-xi(2))*(1.+xi(3)))/8.; dNdxi(6,1) = ((1.-xi(2))*(1.+xi(3))*(xi(1)-xi(2)+xi(3)-2.)+(1.+xi(1))*(1.-xi(2))*(1.+xi(3)))/8.; dNdxi(6,2) = (-(1.+xi(1))*(1.+xi(3))*(xi(1)-xi(2)+xi(3)-2.)-(1.+xi(1))*(1.-xi(2))*(1.+xi(3)))/8.; dNdxi(6,3) = ((1.+xi(1))*(1.-xi(2))*(xi(1)-xi(2)+xi(3)-2.)+(1.+xi(1))*(1.-xi(2))*(1.+xi(3)))/8.; dNdxi(7,1) = ((1.+xi(2))*(1.+xi(3))*(xi(1)+xi(2)+xi(3)-2.)+(1.+xi(1))*(1.+xi(2))*(1.+xi(3)))/8.; dNdxi(7,2) = ((1.+xi(1))*(1.+xi(3))*(xi(1)+xi(2)+xi(3)-2.)+(1.+xi(1))*(1.+xi(2))*(1.+xi(3)))/8.; dNdxi(7,3) = ((1.+xi(1))*(1.+xi(2))*(xi(1)+xi(2)+xi(3)-2.)+(1.+xi(1))*(1.+xi(2))*(1.+xi(3)))/8.; dNdxi(8,1) = (-(1.+xi(2))*(1.+xi(3))*(-xi(1)+xi(2)+xi(3)-2.)-(1.-xi(1))*(1.+xi(2))*(1.+xi(3)))/8.; dNdxi(8,2) = ((1.-xi(1))*(1.+xi(3))*(-xi(1)+xi(2)+xi(3)-2.)+(1.-xi(1))*(1.+xi(2))*(1.+xi(3)))/8.; dNdxi(8,3) = ((1.-xi(1))*(1.+xi(2))*(-xi(1)+xi(2)+xi(3)-2.)+(1.-xi(1))*(1.+xi(2))*(1.+xi(3)))/8.; dNdxi(9,1) = -2.*xi(1)*(1.-xi(2))*(1.-xi(3))/4.; dNdxi(9,2) = -(1.-xi(1)^2)*(1.-xi(3))/4.; dNdxi(9,3) = -(1.-xi(1)^2)*(1.-xi(2))/4.; dNdxi(10,1) = (1.-xi(2)^2)*(1.-xi(3))/4.; dNdxi(10,2) = -2.*xi(2)*(1.+xi(1))*(1.-xi(3))/4.; dNdxi(10,3) = -(1.-xi(2)^2)*(1.+xi(1))/4.; dNdxi(11,1) = -2.*xi(1)*(1.+xi(2))*(1.-xi(3))/4.; dNdxi(11,2) = (1.-xi(1)^2)*(1.-xi(3))/4.; dNdxi(11,3) = -(1.-xi(1)^2)*(1.+xi(2))/4.; dNdxi(12,1) = -(1.-xi(2)^2)*(1.-xi(3))/4.; dNdxi(12,2) = -2.*xi(2)*(1.-xi(1))*(1.-xi(3))/4.; dNdxi(12,3) = -(1.-xi(2)^2)*(1.-xi(1))/4.; dNdxi(13,1) = -2.*xi(1)*(1.-xi(2))*(1.+xi(3))/4.; dNdxi(13,2) = -(1.-xi(1)^2)*(1.+xi(3))/4.; dNdxi(13,3) = (1.-xi(1)^2)*(1.-xi(2))/4.; dNdxi(14,1) = (1.-xi(2)^2)*(1.+xi(3))/4.; dNdxi(14,2) = -2.*xi(2)*(1.+xi(1))*(1.+xi(3))/4.; dNdxi(14,3) = (1.-xi(2)^2)*(1.+xi(1))/4.; dNdxi(15,1) = -2.*xi(1)*(1.+xi(2))*(1.+xi(3))/4.; dNdxi(15,2) = (1.-xi(1)^2)*(1.+xi(3))/4.; dNdxi(15,3) = (1.-xi(1)^2)*(1.+xi(2))/4.; dNdxi(16,1) = -(1.-xi(2)^2)*(1.+xi(3))/4.; dNdxi(16,2) = -2.*xi(2)*(1.-xi(1))*(1.+xi(3))/4.; dNdxi(16,3) = (1.-xi(2)^2)*(1.-xi(1))/4.; dNdxi(17,1) = -(1.-xi(2))*(1.-xi(3)^2)/4.; dNdxi(17,2) = -(1.-xi(1))*(1.-xi(3)^2)/4.; dNdxi(17,3) = -xi(3)*(1.-xi(1))*(1.-xi(2))/2.; dNdxi(18,1) = (1.-xi(2))*(1.-xi(3)^2)/4.; dNdxi(18,2) = -(1.+xi(1))*(1.-xi(3)^2)/4.; dNdxi(18,3) = -xi(3)*(1.+xi(1))*(1.-xi(2))/2.; dNdxi(19,1) = (1.+xi(2))*(1.-xi(3)^2)/4.; dNdxi(19,2) = (1.+xi(1))*(1.-xi(3)^2)/4.; dNdxi(19,3) = -xi(3)*(1.+xi(1))*(1.+xi(2))/2.; dNdxi(20,1) = -(1.+xi(2))*(1.-xi(3)^2)/4.; dNdxi(20,2) = (1.-xi(1))*(1.-xi(3)^2)/4.; dNdxi(20,3) = -xi(3)*(1.-xi(1))*(1.+xi(2))/2.; end end end