a 2D version
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

75 lines
2.2 KiB

4 years ago
%--------------------------
% @Author: Jingqiao Hu
% @Date: 2020-10-29 15:04:18
% @LastEditTime: 2021-01-12 15:09:02
%--------------------------
function [deform1, converged] = line_search(rhs, du, fext, deform0, ref_nodes, ...
dN_c, gradN, dF, edofMat, fixeddofs, lx, u, l, opt_rotation)
maxloop = 10; loop = 0;
nodeNum = size(deform0, 1);
deform0 = reshape(deform0', [], 1);
ref0 = reshape(ref_nodes', [], 1);
step = 1;
approxE0 = -sum(dot(rhs, deform0 - ref0));
norm0 = max(max(abs(rhs)));
found_step = 0;
while loop < maxloop
loop = loop+1;
deform1 = step*du + deform0;
deform2 = reshape(deform1,2,[])';
U = reshape(deform2' - ref_nodes', [], 1);
if opt_rotation
R = update_corotation_matrix(edofMat, U, dN_c);
[F, dF] = corotational_F(R, edofMat, ref_nodes, deform2, gradN);
else
F = deform_grad_1scale(edofMat, U, gradN); % {4}: 4*eleNum
end
% F = deform_gradient(eles, reshape(deform1,2,[])', ref_nodes, gradN); % {4*eleNum}: 2*2
inverted = check_inversion(F);
step = step*0.5;
if (inverted)
continue;
% regard as not converged,
% it is advised that the deformation gradient F be temporarily replaced by
% the nearest physically plausible value ~F
end
fint = elastic_force(F, dF, nodeNum, edofMat, lx, u, l);
rhs = fint + fext;
rhs(fixeddofs) = 0;
approxE1 = -sum(dot(rhs, deform1 - ref0));
norm1 = max(max(abs(rhs)));
if( (norm1<norm0) || (approxE1<approxE0) )
found_step=1;
break;
end
end
converged = 0;
if(~found_step || (abs(norm1-norm0)<1e-5))
converged = 1;
end
deform1(fixeddofs) = deform0(fixeddofs);
deform1 = reshape(deform1, 2, [])';
end
function inverted = check_inversion(F)
inverted = false;
eps = 1e-5;
for gp = 1:4
F1 = F{gp, 1}; % 4 * eleNum
Fdet = F1(1, :).*F1(4, :) - F1(2, :).*F1(3, :);
if (sum(Fdet <= eps) > 0)
inverted = true;
return;
end
end
end