nks221
Student
- Sep 5, 2021
- 1
Main file starts from %1: showing error - Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.
% 1. Remodelling procedure
for i = 1:N_ele
if (emat(i,a)<101)
% CALCULATION OF DELTA RHO AND DELTA MAT
% dr (u,rho,k,B,s,adapt)
% u strain energy density
% rho density
% k reference stimulus for proximal femur
% B remodeling ratio coefficient
% s lazy zone parameter
% adapt = 0 --> step adaptation
% adapt = 1 --> linear adaptation
[d_rho(i,a), d_mat(i,a)] = dr(sed(i,a), emat(i,a)*0.02, 0.004, B(....), 0.35, 1);
% MATERIAL UPDATE
emat(i,a+1) = emat(i,a) + d_mat(i,a);
% MATERIAL CONTROL
if (emat(i,a+1)>100)
emat(i,a+1) = 100;
elseif (emat(i,a+1)<4)
emat(i,a+1) = 4;
end
else
emat(i,a+1) = emat(i,a);
end
end
function file is called in line ([d_rho(i,a), d_mat(i,a)] = dr(sed(i,a), emat(i,a)0.02, 0.004, B(....), 0.35, 1)
which is described below
% 1. Step adaptation
Function [d_rho,d_mat] = dr(u,rho,k,B,s,adapt)
% u strain energy density [MPa]
% rho density in [g/cm3]
% k reference stimulus
% B remodeling coefficiant
% s range of lazy zone
% adapt=0 step adaptation
% adapt=1 linear adaptation
d_rho = 0;
stimu = u/rho; % Stimulus
% Adaptation functions
if (adapt == 0)
% step adaptation
if stimu > k*(1+s)
d_mat = 1;
d_rho = 0.02;
elseif stimu < k*(1-s)
d_mat = -1;
d_roh = -0.02;
else
d_roh = 0;
d_mat = 0;
end
% 2. Linear adaptation.
elseif (adapt == 1)
% linear adaptation
if stimu > k*(1+s)
d_rho = B*(u/rho-k*(1+s));
d_mat = round(d_rho/0.02);
if d_mat > 4
d_mat = 4;
end
elseif stimu < k*(1-s)
d_rho = B*(u/rho-k*(1-s));
d_mat = floor(d_rho/0.02);
if d_mat < -6
d_mat = -6;
end
else
d_rho = 0;
d_mat = 0;
end
else
sprintf('Wrong adaptation!')
stop;
end
% 1. Remodelling procedure
for i = 1:N_ele
if (emat(i,a)<101)
% CALCULATION OF DELTA RHO AND DELTA MAT
% dr (u,rho,k,B,s,adapt)
% u strain energy density
% rho density
% k reference stimulus for proximal femur
% B remodeling ratio coefficient
% s lazy zone parameter
% adapt = 0 --> step adaptation
% adapt = 1 --> linear adaptation
[d_rho(i,a), d_mat(i,a)] = dr(sed(i,a), emat(i,a)*0.02, 0.004, B(....), 0.35, 1);
% MATERIAL UPDATE
emat(i,a+1) = emat(i,a) + d_mat(i,a);
% MATERIAL CONTROL
if (emat(i,a+1)>100)
emat(i,a+1) = 100;
elseif (emat(i,a+1)<4)
emat(i,a+1) = 4;
end
else
emat(i,a+1) = emat(i,a);
end
end
function file is called in line ([d_rho(i,a), d_mat(i,a)] = dr(sed(i,a), emat(i,a)0.02, 0.004, B(....), 0.35, 1)
which is described below% 1. Step adaptation
Function [d_rho,d_mat] = dr(u,rho,k,B,s,adapt)
% u strain energy density [MPa]
% rho density in [g/cm3]
% k reference stimulus
% B remodeling coefficiant
% s range of lazy zone
% adapt=0 step adaptation
% adapt=1 linear adaptation
d_rho = 0;
stimu = u/rho; % Stimulus
% Adaptation functions
if (adapt == 0)
% step adaptation
if stimu > k*(1+s)
d_mat = 1;
d_rho = 0.02;
elseif stimu < k*(1-s)
d_mat = -1;
d_roh = -0.02;
else
d_roh = 0;
d_mat = 0;
end
% 2. Linear adaptation.
elseif (adapt == 1)
% linear adaptation
if stimu > k*(1+s)
d_rho = B*(u/rho-k*(1+s));
d_mat = round(d_rho/0.02);
if d_mat > 4
d_mat = 4;
end
elseif stimu < k*(1-s)
d_rho = B*(u/rho-k*(1-s));
d_mat = floor(d_rho/0.02);
if d_mat < -6
d_mat = -6;
end
else
d_rho = 0;
d_mat = 0;
end
else
sprintf('Wrong adaptation!')
stop;
end
