Homec4science
Diffusion MSM (master) Browse

msm/ShoulderModel_HHT/postprocessing_contact_pattern_HH.m
73754adfee53master

Last Change
Raw File

postprocessing_contact_pattern_HH.m
View Options

function postprocessing_contact_pattern_HH(CMDATA, PPDATA, IDDATA)
%{
Function for defining the contact pattern on the humeral head.
--------------------------------------------------------------------------
Syntax :
postprocessing_contact_pattern_HH(CMDATA, PPDATA, IDDATA)
--------------------------------------------------------------------------
File Description :
This function plots the contact pattern on the humeral head. To this end, we
wrap the results defined for the amount of penetration on the humeral head
surface.
--------------------------------------------------------------------------
%}
theta = [0:CMDATA.SpanStepSize:2*pi]; % longtitude 0<theta<2pi on the xy from x
phi = [0:CMDATA.SpanStepSize:pi]; % latitude 0<phi<pi from z
height = [0:5e-3:100e-3]; % the height of the humerus (for illustration only)
% points to draw the sphere
% ------------------------------------------------------------------------
% define the top and bottom points of the sphere (and their associated
% surafce areas)
X(1,1) = CMDATA.r*cos(theta(1))*sin(phi(1));
Y(1,1) = CMDATA.r*cos(phi(1));
Z(1,1) = CMDATA.r*sin(theta(1))*sin(phi(1));
X(size(theta,2),size(phi,2)) = CMDATA.r*cos(theta(1))*sin(phi(size(phi,2)));
Y(size(theta,2),size(phi,2)) = CMDATA.r*cos(phi(size(phi,2)));
Z(size(theta,2),size(phi,2)) = CMDATA.r*sin(theta(1))*sin(phi(size(phi,2)));
% span the whole surface of the humeral head and define the center points
% of the mesh grids. And rightafter their associated surace areas.
for i = 1:size(theta,2)-1 % theta = 0 and 2pi result in the same points
for j = 2:size(phi,2)-1 % exclude the top and bottom points
X(i,j) = CMDATA.r*cos(theta(i))*sin(phi(j));
Y(i,j) = CMDATA.r*cos(phi(j));
Z(i,j) = CMDATA.r*sin(theta(i))*sin(phi(j));
end
end
% points to draw the cylinder
% ------------------------------------------------------------------------
for i = 1:size(theta,2)-1
c1_b(1,i) = CMDATA.r_h*cos(theta(i));
c3_b(1,i) = CMDATA.r_h*sin(theta(i));
end
% ------------------------------------------------------------------------
figure % the main figure of this function
for TimeStep = 1:size(IDDATA.t,2)-1 % to rotate the arm the rotation matrix is considered
% at each time step for the whole time
% horizon
u = [IDDATA.Coordinate_Evolution(:,TimeStep);
IDDATA.Coordinate_Evolution_D(:,TimeStep)]; % the generalized coordinate and velocity vector
% u = [0;IDDATA.ImposedMotion(4,TimeStep);0;0;0;0]; % the generalized coordinate and velocity vector
Rot_H2T = versatile_function('define Rot_H2T',u); % define the rotation matrix
% rotate the sphere---------------------------------------------------
% define the top and bottom points of the sphere (and their associated
% surafce areas)
Xr(1,1) = Rot_H2T(1,:)*[X(1,1) Y(1,1) Z(1,1)]'; % rotated H
Yr(1,1) = Rot_H2T(2,:)*[X(1,1) Y(1,1) Z(1,1)]';
Zr(1,1) = Rot_H2T(3,:)*[X(1,1) Y(1,1) Z(1,1)]';
Xr(size(theta,2),size(phi,2)) = Rot_H2T(1,:)*[X(size(theta,2),size(phi,2)) Y(size(theta,2),size(phi,2)) Z(size(theta,2),size(phi,2))]'; % rotated H
Yr(size(theta,2),size(phi,2)) = Rot_H2T(2,:)*[X(size(theta,2),size(phi,2)) Y(size(theta,2),size(phi,2)) Z(size(theta,2),size(phi,2))]';
Zr(size(theta,2),size(phi,2)) = Rot_H2T(3,:)*[X(size(theta,2),size(phi,2)) Y(size(theta,2),size(phi,2)) Z(size(theta,2),size(phi,2))]';
% span the whole surface of the humeral head and define the center points
% of the mesh grids. And rightafter their associated surace areas.
for i = 1:size(theta,2)-1 % theta = 0 and 2pi result in the same points
for j = 2:size(phi,2)-1 % exclude the top and bottom points
Xr(i,j) = Rot_H2T(1,:)*[X(i,j) Y(i,j) Z(i,j)]'; % rotated H
Yr(i,j) = Rot_H2T(2,:)*[X(i,j) Y(i,j) Z(i,j)]';
Zr(i,j) = Rot_H2T(3,:)*[X(i,j) Y(i,j) Z(i,j)]';
end
end
% rotate the cylinder-------------------------------------------------
for j = 1:size(height,2)
for i= 1:size(theta,2)-1
c1_br(j,i) = Rot_H2T(1,:)*[c1_b(1,i) -height(j) c3_b(1,i)]';
c3_br(j,i) = Rot_H2T(3,:)*[c1_b(1,i) -height(j) c3_b(1,i)]';
c2_br(j,i) = Rot_H2T(2,:)*[c1_b(1,i) -height(j) c3_b(1,i)]';
end
end
surf(Xr,Zr,Yr,1e3*PPDATA.deformation{TimeStep}) % plot the sphere together
% with the deformations magnitudes
% view(-120,5) % http://ch.mathworks.com/help/matlab/ref/view.html
hold on
surf(c1_br,c3_br,c2_br,zeros(size(c1_br,1),size(c1_br,2))) % plot the cylinder
cmap = jet(50); % define the color map of the figures
cmap(1,:)=[0.6000 0.6000 0.6000]; % define the base color of the bones
colormap(cmap);
h1 = colorbar;
set(h1,'fontsize',16); % set the font of the colorbar
ylabel(h1, 'cartilage deformation [mm]')
% show extra information about the contact
str1 = ['abduction angle = ',num2str(-(180/pi)*IDDATA.ImposedMotion(4,TimeStep)), ' [deg]']; % show the abduction angle on the plot
str1_prop = annotation('textbox',[0.01 0.01 0.5 0.05],'String',str1,'FitBoxToText','off','BackgroundColor','w');
str1_prop.FontSize = 16;
str2 = ['contact area = ',num2str(PPDATA.cont_area(TimeStep)), ' [mm\^{2}]']; % show the contact area on the plot
str2_prop = annotation('textbox',[0.01 0.07 0.5 0.05],'String',str2,'FitBoxToText','off','BackgroundColor','w');
str2_prop.FontSize = 16;
title('contact pattern on the humeral head',...
'fontsize', 11, ...
'fontweight', 'bold', ...
'fontname', 'sansserif');
% make the animation
pause(0.3)
if TimeStep < size(IDDATA.t,2)-1
cla
else
% keep the last figure
end
% whitebg([1,1,1]);
axis equal
axis off
end
return;
c4science · Help