Homec4science
Diffusion MSM (command_line) Browse

msm/ShoulderModel_HHT/postprocessing_contact_model.m
73754adfee53command_line

Last Change
Raw File

postprocessing_contact_model.m
View Options

function PPDATA = postprocessing_contact_model(TimeStep,Rot_H2T,DRot_H2T,r_HH,Dr_HH,CMDATA, PPDATA)
%{
Function for defining the contacts data at each joint configuration.
--------------------------------------------------------------------------
Syntax :
PPDATA = postprocessing_contact_model(TimeStep,Rot_H2T,DRot_H2T,r_HH,Dr_HH,CMDATA, PPDATA)
--------------------------------------------------------------------------
Function Description :
This function is responsible for postprocessing the achieved results and
contructs the PPDATA structure.
The function receives the current joint configuration (rotation matrix and
the humeral head translation vector as well as their associated rate of
changes) and defines the contact force at each point on the humeral head
surface based on whether the point is in contact or not. We already defined
the points in the humerus body-fixed coordinate frame through
initilize_contact_model(.) saved in CMDATA.
This function is a replication of the "contact_model" to reproduce the
contact data based on the resulted generalized coordinates and velocities.
Here are the outputs of the function:
1- deformation: based on the points of the humeral head.
2- F_contact: contact forces at each points of the humeral
head.
3- Dir_contact: directions of the contact forces at each
points of the humeral head.
4- F_resultant: the resultant force of all the forces acting on
each of the points on the humeral head. A single vector for
each time step.
5- F_resultant_magnitude: the magnitude of the resultant force
at each time step. The magnitude of the force achieved form 4.
6- h_inT: the positions of the points on the humeral head for
which the penetration is not zero. These positions are defined
in the thorax coordinate frame. Those points that are not in
contact were set to zero.
7- cont_area: the contact area of the points that are in
contact. The other values are set to zero.
8- max_deformation: the maximum cartilage deformation at each
time step.
--------------------------------------------------------------------------
%}
h = CMDATA.h; % receive h_i s in the frame {H}
a = CMDATA.a; % receive associtaed surface area of the h_i s
for i = 1:size(h,1)-1 % associated with longtitude (theta)
for j = 2:size(h,2)-1 % associated with latitude (phi)
h_inT{i,j} = r_HH + Rot_H2T*h{i,j}; % hi in the {T} frame
u(i,j) = norm(h_inT{i,j},2) - CMDATA.r_g; % the penetration of the point hi
% in the glenoid fossa
Dh_inT{i,j} = Dr_HH + DRot_H2T*h{i,j}; % rate of change of hi
Du(i,j) = dot(h_inT{i,j},Dh_inT{i,j})/norm(h_inT{i,j},2); % rate of change of the ui
if gt(u(i,j),0) & lt(u(i,j),0.98*CMDATA.b) % contact condition of the points
F_contact(i,j) = -a(i,j)*CMDATA.s*log(1-u(i,j)/CMDATA.b) + CMDATA.c*Du(i,j);
Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
Rho_contact{i,j} = Rot_H2T*h{i,j};
h_inT{i,j} = h_inT{i,j};
cont_area(i,j) = a(i,j);
elseif ge(u(i,j),0.98*CMDATA.b)
%u(i,j) = 0.98*b;
F_contact(i,j) = -a(i,j)*CMDATA.s*log(1-0.98) + CMDATA.c*Du(i,j);
Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
Rho_contact{i,j} = Rot_H2T*h{i,j};
u(i,j) = 0.98*CMDATA.b;
h_inT{i,j} = h_inT{i,j};
cont_area(i,j) = a(i,j);
%disp(['u' num2str(i) ',' num2str(i) '=' num2str(u(i,j)) '>=0.98b']);
else % and if there is no contact there will be no contact force
F_contact(i,j) = 0;
Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
Rho_contact{i,j} = Rot_H2T*h{i,j};
u(i,j) = 0;
h_inT{i,j} = 0;
cont_area(i,j) = 0;
end
end
end
% the same calculations for the top and bottom points
% THE TOP POINTS
h_inT{1,1} = r_HH + Rot_H2T*h{1,1}; % hi in the {T} frame
u(1,1) = norm(h_inT{1,1},2) - CMDATA.r_g; % the penetration of the point hi
% in the glenoid fossa
Dh_inT{1,1} = Dr_HH + DRot_H2T*h{1,1}; % rate of change of h1
Du(1,1) = dot(h_inT{1,1},Dh_inT{1,1})/norm(h_inT{1,1},2); % rate of change of the ui
if gt(u(1,1),0) & lt(u(1,1),0.98*CMDATA.b) % contact condition of the points
F_contact(1,1) = -a(1,1)*CMDATA.s*log(1-u(1,1)/CMDATA.b) + CMDATA.c*Du(1,1);
Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
Rho_contact{1,1} = Rot_H2T*h{1,1};
h_inT{1,1} = h_inT{1,1};
cont_area(1,1) = a(1,1);
elseif ge(u(1,1),0.98*CMDATA.b)
%u(i,j) = 0.98*b;
F_contact(1,1) = -a(1,1)*CMDATA.s*log(1-0.98) + CMDATA.c*Du(1,1);
Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
Rho_contact{1,1} = Rot_H2T*h{1,1};
u(1,1) = 0.98*CMDATA.b;
h_inT{1,1} = h_inT{1,1};
cont_area(1,1) = a(1,1);
%disp(['u' num2str(1) ',' num2str(1) '=' num2str(u(1,1)) '>=0.98b']);
else % and if there is no contact there will be no contact force
F_contact(1,1) = 0;
Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
Rho_contact{1,1} = Rot_H2T*h{1,1};
u(1,1) = 0;
h_inT{1,1} = 0;
cont_area(1,1) = 0;
end
% THE BOTTOM POINT
k = size(h,1); l = size(h,2);
h_inT{k,l} = r_HH + Rot_H2T*h{k,l}; % hi in the {T} frame
u(k,l) = norm(h_inT{k,l},2) - CMDATA.r_g; % the penetration of the point hi
% in the glenoid fossa
Dh_inT{k,l} = Dr_HH + DRot_H2T*h{k,l}; % rate of change of h bottom
Du(k,l) = dot(h_inT{k,l},Dh_inT{k,l})/norm(h_inT{k,l},2); % rate of change of the ui
if gt(u(k,l),0) & lt(u(k,l),0.98*CMDATA.b) % contact condition of the points
F_contact(k,l) = -a(k,l)*CMDATA.s*log(1-u(k,l)/CMDATA.b) + CMDATA.c*Du(k,l);
Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
Rho_contact{k,l} = Rot_H2T*h{k,l};
h_inT{k,l} = h_inT{k,l};
cont_area(k,l) = a(k,l);
elseif ge(u(k,l),0.98*CMDATA.b)
%u(k,l) = = 0.98*b;
F_contact(k,l) = -a(k,l)*CMDATA.s*log(1-0.98) + CMDATA.c*Du(k,l);
Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
Rho_contact{k,l} = Rot_H2T*h{k,l};
u(k,l) = 0.98*CMDATA.b;
h_inT{k,l} = h_inT{k,l};
cont_area(k,l) = a(k,l);
%disp(['u' num2str(k) ',' num2str(l) '=' num2str(u(k,l)) '>=0.98b']);
else % and if there is no contact there will be no contact force
F_contact(k,l) = 0;
Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
Rho_contact{k,l} = Rot_H2T*h{k,l};
u(k,l) = 0;
h_inT{k,l} = 0;
cont_area(k,l) = 0;
end
% NOTE: these values are defined for the inputted time step. We make the structure
% containing the information for the whole time horizon in the for loop in
% the main function of the postprocessing.
PPDATA.deformation{TimeStep} = u; % deformations of all the points on the humeral head
PPDATA.max_deformation(TimeStep) = max(u(:)); % maximum deformations at each time step
PPDATA.F_contact{TimeStep} = F_contact;
PPDATA.Dir_contact{TimeStep} = Dir_contact;
PPDATA.h_inT{TimeStep} = h_inT;
PPDATA.cont_area(TimeStep) = 1e6*sum(cont_area(:)); % contact area [mm2]
% % define the resultant contact force
F_resultant = postprocessing_contact_force(Dir_contact,F_contact); % function to deine the resultant contact force
PPDATA.F_resultant{TimeStep} = F_resultant;
PPDATA.F_resultant_magnitude(TimeStep) = norm(F_resultant,2);% magnitude of the resultant muscle force
% deformation = u; % deformations of all the points on the humeral head
% F_resultant = postprocessing_contact_force(Dir_contact,F_contact); % function to deine the resultant contact force
% F_resultant_magnitude = norm(F_resultant,2); % magnitude of the resultant muscle force
return;
c4science · Help