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

addpath('..\two_dim_fea');

width = 30; height = 30; h = 1;
m = 10; n = 10;
hw = width/m; hh = height/n;

E = 10.0; nu = 0.48;
c = 0.04*E;
rho = 1.00;
Epfloor = 0.02; cpfloor = 0;
list = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ];

tp = 0.5e+3; %tp = 1e+3;
vp = 0.8*(height/3)/tp;
th = 0.5e+3; %th = 1e+3;
tf = 0.5e+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,2*m);
for k=1:2*m; A_fixed(k,k)=1; end

% 押す節点
A_push = zeros(npoint*2,8); % P94 - P97
for k=1:8; A_push(94*2-1+(k-1),k) = 1; end

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

% 上部を押している
square_object_push = @(t,q) square_object_jump_push_param(t,q,npoint,M,B,K,4,A_push,push,dotpush,m,A_fixed,alpha);
interval = [0, tp+th];
qinit = zeros(npoint*4,1);
[time_push, q_push] = ode23tb(square_object_push, interval, qinit);

% 上部を解放した後
external_force = @(t,n,un,vn) floor_force_param(t,npoint,un,vn, list,h,hw,Epfloor,cpfloor); % 下端の床反力
square_object_free = @(t,q) square_object_jump_free_param(t,q,npoint,M,B,K,external_force);
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);
%figure('Position', [403 246 560 420]);
figure('Position', [403 246 210 400]);

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,90]);
    xticks([-10:10:40])
    yticks([-10:10:90])
    pbaspect([1 2 1]);
    grid on;
    filename = strcat('10-10/deform-', num2str(floor(t),'%04d'), '.png');
    saveas(gcf, filename, 'png');
end

figure('Position', [403 246 560 420]);
clf('reset');
ts = time(1);
te = time(end);
fr = 1;
clear M;
for t = 0:0.1*100:tp+th+tf
    index = nearest_index(time, t);
    draw_triangles(point, triangle, q(index,1:npoint*2));
    hold off;
    xlim([-10,40]);
    ylim([-10,90]);
    xticks([-10:10:40])
    yticks([-10:10:90])
    pbaspect([1 2 1]);
    title(['time ' num2str(t/1000,"%3.2f")]);
    grid on;
    drawnow;
    M(fr) = getframe(gcf);
    fr = fr + 1;
    disp(t/1000);
end
M(fr) = getframe(gcf);

v = VideoWriter('two_dim_jump_10_10');
open(v);
writeVideo(v, M);
close(v);
