Skip to content
Snippets Groups Projects
Commit f5c00b32 authored by Boman Romain's avatar Boman Romain
Browse files

add more comments

parent e13e7082
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
nastran/Boneload_version_6_w_4_methods_ZJT.m
+ 43
− 29
View file @ f5c00b32
......@@ -29,8 +29,8 @@
% Focal_nodes: focal point coordinates
% for each group size = (nbgroups,3)
%
% MG_{i}_Focal_Point: shortcut to 'Focal_nodes' for group i
% MG_{i}_Muscle_Force: shortcut to 'Force' for group i
% MG_{i}_Focal_Point: shortcut to 'Focal_nodes' for group i size=(3)
% MG_{i}_Muscle_Force: shortcut to 'Forces' for group i size=(1)
%
% method: loading method
% 'uniform_traction': "ad hoc uniform traction model"
......@@ -50,8 +50,8 @@
%
% MG_{i}_Element_Area size = (nbelem_i)
% MG_{i}_Elem_Centroid size = (nbelem_i, 3)
% MG_{i}_Traction size = (nbelem_i, 3)?
% MG_{i}_Normal size = (nbelem_i, 3)
% MG_{i}_Traction : tangential comp of the force size = (nbelem_i, 3)
% MG_{i}_Normal : normal comp of the force size = (nbelem_i, 3)
% MG_{i}_free_elem_edges
% MG_{i}_outter_edge
% MG_{i}_Outter_Edge_Nodes: list of boundardy edges size = (nbonod)
......@@ -62,7 +62,11 @@
% Rem_Neg_Mom: flags (0 or 1) used for removing negative moments
% size = (nbgroups)
%
%
% the script also copies the date in a structure while computing the forces
% Muscle_Info.Muscle_attachment_X
% Muscle_Info.Muscle_Traction
% Muscle_Info.Muscle_Elem_Centroid
% ...
%
%
%
......@@ -182,6 +186,7 @@ end
for jj = 1:1:max_number_of_muscle_groups
eval(sprintf('[MG_%d_Element_Area, MG_%d_Elem_Centroid, MG_%d_Traction, MG_%d_Normal] = Element_info(MG_%d_Muscle_attachment_X,MG_%d_Muscle_attachment_Y,MG_%d_Muscle_attachment_Z,MG_%d_Muscle_Force,MG_%d_Focal_Point);',jj,jj,jj,jj,jj,jj,jj,jj,jj))
% [RB] ! MG_%d_Traction => "traction" = tangential component of the force!!
clf
figure(21)
set(gcf,'Position',[500 500 500 100]);
......@@ -343,6 +348,9 @@ for jj = 1:1:max_number_of_muscle_groups
MG_Elem_Normal = Muscle_Info.Muscle_Normal(elem_num,:);
% THIS CHANGES THE PRESSURE MAGNITUDE AND DIRECTION
%
% this is the "mean pressure" - same for all elements
% <----------------------------------------------------------->
Pressure_load(elem_num,:) = (Muscle_Info.Muscle_Force./sum(Muscle_Info.Muscle_Elem_Area))*MG_Elem_CENT_to_FP_unit;
%%% Remove Negative Moments
......@@ -352,13 +360,24 @@ for jj = 1:1:max_number_of_muscle_groups
if dot(Moment_check, [L_TMJ - R_TMJ]./norm([L_TMJ - R_TMJ]))<0
Pressure_load(elem_num,:) = 0*(MG_Elem_CENT_to_FP_unit);
end
end
end
end
% [RB]
% The method consists in computing a (scalar) pressure_mean = F_prescribed/total_Area
% Then the "pressure_load" (3 comp vector with unit=pressure) equal
% pressure_load = pressure_mean*(unit vector from centroid to focal point)
% [RB]: no "normalize_to_input_force" here!
%% Gradient Traction (Pressure) Method
elseif strcmp(method,'gradient_traction') ==1
MG_Elem_CENT_to_FP_dist = zeros(Muscle_Info.num_of_plates,1);
% [RB] same as above except that the force amplitude is scaled by a factor
% that depends on the distance of the facet to the focal point.
% (1 for the farthest, less for nearer elements)
MG_Elem_CENT_to_FP_dist = zeros(Muscle_Info.num_of_plates,1);
for elem_num = 1:1:Muscle_Info.num_of_plates
MG_Elem_CENT_to_FP_dist(elem_num) = norm(Muscle_Info.Focal_Point - Muscle_Info.Muscle_Elem_Centroid (elem_num,:));
......@@ -373,30 +392,21 @@ for jj = 1:1:max_number_of_muscle_groups
MG_Elem_Normal = Muscle_Info.Muscle_Normal(elem_num,:);
% THIS CHANGES THE PRESSURE MAGNITUDE AND DIRECTION
Pressure_load(elem_num,:) = gradient_scale(elem_num).*(Muscle_Info.Muscle_Force./sum(Muscle_Info.Muscle_Elem_Area))*MG_Elem_CENT_to_FP_unit;
% %% PROGRESS BAR
% if mod(elem_num,100)==0
% elem_num
% % computation here %
% waitbar(elem_num/Muscle_Info.num_of_plates,h_bar)
% end
Pressure_load(elem_num,:) = gradient_scale(elem_num).*(Muscle_Info.Muscle_Force./sum(Muscle_Info.Muscle_Elem_Area))*MG_Elem_CENT_to_FP_unit;
end
normalize_to_input_force = Muscle_Info.Muscle_Force./sum(sqrt(Pressure_load(:,1).^2+Pressure_load(:,2).^2+Pressure_load(:,3).^2).*(Muscle_Info.Muscle_Elem_Area(:)));
normalize_to_input_force = Muscle_Info.Muscle_Force./sum(sqrt(Pressure_load(:,1).^2+Pressure_load(:,2).^2+Pressure_load(:,3).^2).*(Muscle_Info.Muscle_Elem_Area(:)));
for elem_num = 1:1:Muscle_Info.num_of_plates
%%% Remove Negative Moments
% REMOVE THE PRESSURE THAT CAUSES A NEGATIVE MOMENT ABOUT THE TMJ
if Rem_Neg_Mom(jj) ==1
Moment_check = cross(Muscle_Info.Muscle_Elem_Centroid (elem_num,:) - L_TMJ,Pressure_load(elem_num,:));
if dot(Moment_check, [L_TMJ - R_TMJ]./norm([L_TMJ - R_TMJ]))<0
Pressure_load(elem_num,:) = 0*(MG_Elem_CENT_to_FP_unit);
end
for elem_num = 1:1:Muscle_Info.num_of_plates
%%% Remove Negative Moments
% REMOVE THE PRESSURE THAT CAUSES A NEGATIVE MOMENT ABOUT THE TMJ
if Rem_Neg_Mom(jj) ==1
Moment_check = cross(Muscle_Info.Muscle_Elem_Centroid (elem_num,:) - L_TMJ,Pressure_load(elem_num,:));
if dot(Moment_check, [L_TMJ - R_TMJ]./norm([L_TMJ - R_TMJ]))<0
Pressure_load(elem_num,:) = 0*(MG_Elem_CENT_to_FP_unit);
end
end
end
%% TANGENTIAL ONLY
......@@ -411,6 +421,7 @@ for jj = 1:1:max_number_of_muscle_groups
if dot(MG_Elem_CENT_to_FP_unit,MG_Elem_Normal) > 0
Pressure_load(elem_num,:) = (Muscle_Info.Muscle_Force./sum(Muscle_Info.Muscle_Elem_Area))*MG_Elem_CENT_to_FP_unit;
else
% tv = "tangent vector"
tv = Muscle_Info.Muscle_Traction(elem_num,:)./(norm(Muscle_Info.Muscle_Traction(elem_num,:)));
Pressure_load(elem_num,:) = (Muscle_Info.Muscle_Force./sum(Muscle_Info.Muscle_Elem_Area))*(tv);
end
......@@ -423,8 +434,9 @@ for jj = 1:1:max_number_of_muscle_groups
Pressure_load(elem_num,:) = 0*(MG_Elem_CENT_to_FP_unit);
end
end
end
% [RB]: no "normalize_to_input_force" here!
elseif strcmp(method,'tangential_plus_normal_traction') ==1
......@@ -458,7 +470,7 @@ for jj = 1:1:max_number_of_muscle_groups
% We use this formula:
% Tt (tv1 i + tv2 j + tv3 k) + Tn(nv1 i + nv2 j + nv3 k);
% where:
% t1,t2,t3 indicate the direction of the inplane tration vector
% t1,t2,t3 indicate the direction of the inplane traction vector ([RB]: "traction" means "tangent" here)
% n1,n2,n3 indicate the direction of the outward normal vector
% Tt = The Tangential Traction (F/A)
% Tn = Tt*(s(r)/R(r))
......@@ -468,7 +480,9 @@ for jj = 1:1:max_number_of_muscle_groups
% 1069-1088
% FORCE VECTOR DIRECTIONS
% tangent vector
tv = Muscle_Info.Muscle_Traction(elem_num,:)./(norm(Muscle_Info.Muscle_Traction(elem_num,:)));
% normal vector
nv = - Muscle_Info.Muscle_Normal(elem_num,:);
% THIS CHANGES THE PRESSURE MAGNITUDE AND DIRECTION
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment