% 正方形物体の動的な弾性変形 (4&times;4)
% g, cm, msec

width = 30; height = 30; h = 1;
m = 4; n = 4;

E = 10.0; nu = 0.48;
c = 0.04*E;
rho = 1.00;

tp = 1e+3;
vp = 0.8*(height/3)/tp;
th = 1e+3;
tf = 3e+3;

alpha = 1e+6;

npoint = m*n;
[point,triangle] = rectangular_object(m,n,width,height);
         
[lambda,mu] = Lame_constants(E, nu);
[lambdav,muv] = Lame_constants(c, nu);

[Jlambda, Jmu] = connection_matrices(point,triangle,h);
K = lambda*Jlambda + mu*Jmu;
B = lambdav*Jlambda + muv*Jmu;
M = inertia_matrix(point,triangle,h,rho);

% 床に固定する節点
A_fixed = zeros(npoint*2,8);
for k=1:8; A_fixed(k,k)=1; end

% 押す節点
A_push = zeros(npoint*2,4);
A_push(27,1) = 1; A_push(28,2) = 1; A_push(29,3) = 1; A_push(30,4) = 1;

push = @(t) push_param (2, t, tp,vp,th);
dotpush = @(t) dotpush_param (2, t, tp,vp,th);

square_object_push = @(t,q) square_object_push_param(t,q,npoint,M,B,K,2,A_push,push,dotpush,4,A_fixed,alpha);
square_object_free = @(t,q) square_object_free_param(t,q,npoint,M,B,K,4,A_fixed,alpha);

% 上部を押している
interval = [0, tp+th];
qinit = zeros(npoint*4,1);
[time_push, q_push] = ode23tb(square_object_push, interval, qinit);

% 上部を解放した後
interval = [tp+th, tp+th+tf];
sz = size(q_push); tn = sz(1);
qinit = q_push(tn,:);
[time_free, q_free] = ode23tb(square_object_free, interval, qinit);

time = [time_push; time_free];
q = [q_push; q_free];

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

for t = 0:0.1*1000:tp+th+tf
    fprintf("時刻 %f\n", t);
    index = nearest_index(time, t);
    draw_triangles(point, triangle, q(index,1:npoint*2));
    hold off;
    xlim([-10,40]);
    ylim([-10,40]);
    pbaspect([1 1 1]);
    grid on;
    filename = strcat('4-4/deform-', num2str(t), '.eps');
    saveas(gcf, filename, 'eps');
end
