Homec4science
Diffusion MSM (master) Browse

msm/ShoulderModel_HHT/contact_model_dominique.m
73754adfee53master

Last Change
Raw File

contact_model_dominique.m
View Options

function [Dir_contact,Rho_contact,F_contact]= contact_model_dominique(Rot_H2T,DRot_H2T,r_HH,Dr_HH,CMDATA)
%{
Function for defining the contact forces at each joint configuration.
--------------------------------------------------------------------------
Syntax :
[Dir_contact,Rho_contact,F_contact]= contact_model(Rot_H2T,r_HH,CMDATA)
--------------------------------------------------------------------------
Function Description :
The function receives the current joint configuration (rotation matrix and
the humeral head translation vector as well as their associated rate of
change) 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().
--------------------------------------------------------------------------
%}
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
I1(i,j) = (2*CMDATA.b/(CMDATA.b-u(i,j))) + (CMDATA.b-u(i,j))^2/CMDATA.b^2;
F_contact(i,j) = a(i,j)*(CMDATA.b/(CMDATA.b-u(i,j)))*(CMDATA.alpha*CMDATA.beta*(2*exp(CMDATA.beta*(I1(i,j) - 3)) - I1(i,j)) + CMDATA.alpha*CMDATA.beta*(CMDATA.b-u(i,j))^2/CMDATA.b^2 - CMDATA.eta*(I1(i,j) - 3)*2*(CMDATA.b - u(i,j))*Du(i,j)/(CMDATA.b^2));
Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
Rho_contact{i,j} = Rot_H2T*h{i,j};
elseif ge(u(i,j),0.98*CMDATA.b)
u(i,j) = 0.98*CMDATA.b;
I1(i,j) = (2*CMDATA.b/(CMDATA.b-u(i,j))) + (CMDATA.b-u(i,j))^2/CMDATA.b^2;
F_contact(i,j) = a(i,j)*(CMDATA.b/(CMDATA.b-u(i,j)))*(CMDATA.alpha*CMDATA.beta*(2*exp(CMDATA.beta*(I1(i,j) - 3)) - I1(i,j)) + CMDATA.alpha*CMDATA.beta*(CMDATA.b-u(i,j))^2/CMDATA.b^2 - CMDATA.eta*(I1(i,j) - 3)*2*(CMDATA.b - u(i,j))*Du(i,j)/(CMDATA.b^2));
Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
Rho_contact{i,j} = Rot_H2T*h{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};
end
end
end
% the same calculations for the top and bottom points
% the top point
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
I1(1,1) = (2*CMDATA.b/(CMDATA.b-u(1,1))) + (CMDATA.b-u(1,1))^2/CMDATA.b^2;
F_contact(1,1) = a(1,1)*(CMDATA.b/(CMDATA.b-u(1,1)))*(CMDATA.alpha*CMDATA.beta*(2*exp(CMDATA.beta*(I1(1,1) - 3)) - I1(1,1)) + CMDATA.alpha*CMDATA.beta*(CMDATA.b-u(1,1))^2/CMDATA.b^2 - CMDATA.eta*(I1(1,1) - 3)*2*(CMDATA.b - u(1,1))*Du(1,1)/(CMDATA.b^2));
Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
Rho_contact{1,1} = Rot_H2T*h{1,1};
elseif ge(u(1,1),0.98*CMDATA.b)
u(1,1) = 0.98*CMDATA.b;
I1(1,1) = (2*CMDATA.b/(CMDATA.b-u(1,1))) + (CMDATA.b-u(1,1))^2/CMDATA.b^2;
F_contact(1,1) = a(1,1)*(CMDATA.b/(CMDATA.b-u(1,1)))*(CMDATA.alpha*CMDATA.beta*(2*exp(CMDATA.beta*(I1(1,1) - 3)) - I1(1,1)) + CMDATA.alpha*CMDATA.beta*(CMDATA.b-u(1,1))^2/CMDATA.b^2 - CMDATA.eta*(I1(1,1) - 3)*2*(CMDATA.b - u(1,1))*Du(1,1)/(CMDATA.b^2));
Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
Rho_contact{1,1} = Rot_H2T*h{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};
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
I1(k,l) = (2*CMDATA.b/(CMDATA.b-u(k,l))) + (CMDATA.b-u(k,l))^2/CMDATA.b^2;
F_contact(k,l) = a(k,l)*(CMDATA.b/(CMDATA.b-u(k,l)))*(CMDATA.alpha*CMDATA.beta*(2*exp(CMDATA.beta*(I1(k,l) - 3)) - I1(k,l)) + CMDATA.alpha*CMDATA.beta*(CMDATA.b-u(k,l))^2/CMDATA.b^2 - CMDATA.eta*(I1(k,l) - 3)*2*(CMDATA.b - u(k,l))*Du(k,l)/(CMDATA.b^2));
Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
Rho_contact{k,l} = Rot_H2T*h{k,l};
elseif ge(u(k,l),0.98*CMDATA.b)
u(k,l) = 0.98*CMDATA.b;
I1(k,l) = (2*CMDATA.b/(CMDATA.b-u(k,l))) + (CMDATA.b-u(k,l))^2/CMDATA.b^2;
F_contact(k,l) = a(k,l)*(CMDATA.b/(CMDATA.b-u(k,l)))*(CMDATA.alpha*CMDATA.beta*(2*exp(CMDATA.beta*(I1(k,l) - 3)) - I1(k,l)) + CMDATA.alpha*CMDATA.beta*(CMDATA.b-u(k,l))^2/CMDATA.b^2 - CMDATA.eta*(I1(k,l) - 3)*2*(CMDATA.b - u(k,l))*Du(k,l)/(CMDATA.b^2));
Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
Rho_contact{k,l} = Rot_H2T*h{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};
end
%imagesc(u);
%colorbar;
return
c4science · Help