% Comparing calculations based on Cauchy and Green strains
% left end fixed and external force acts downward to the center of right end

% コーシー歪みとグリーン歪みによる計算を比較
% ビームの左端を固定．右端の中央に下向きに力を加える．

% g, cm, sec

addpath('../two_dim_fea');

m = 11; n = 3;
thickness = 1;

[ points, triangles, rectangles ] = rectangular_object( m, n, 10, 2 );

% E = 0.1 MPa; rho = 0.020 g/cm^2
Young = 1.0*1e+6; nu = 0.48; density = 0.020;
% 10.0 N --> fpush = 10.0*1e+5;
force = [0; -1.2*1e+5];

[ lambda, mu ] = Lame_constants( Young, nu );

npoints = size(points,2);
nrectangles = size(rectangles,1);
elastic = Body(npoints, points, [], [], thickness);
elastic = elastic.rectangle_elements(nrectangles, rectangles);
elastic = elastic.mechanical_parameters(density, lambda, mu);

disps_natural = zeros(2,npoints);
un_natural = reshape(disps_natural, [2*npoints,1]);

internal_energy = @(un) internal_energy_params( un, npoints, elastic, force );
internal_energy_Green_strain = @(un) internal_energy_Green_strain_params( un, npoints, elastic, force );

% fix point 1 1+m 1+2m
% a set of constraints: A' un = b (each column specifies one constraint)
A = zeros(2*npoints, 6); b = zeros(6,1);
k = 1;     A(2*k-1, 1) = 1; A(2*k, 2) = 1;
k = 1+m;   A(2*k-1, 3) = 1; A(2*k, 4) = 1;
k = 1+2*m; A(2*k-1, 5) = 1; A(2*k, 6) = 1;

options = optimset('Display','iter', 'MaxFunEvals',3*10^5);
[un_Cauchy, Emin_Cauchy] = fmincon(internal_energy, un_natural, [], [], A', b, [], [], [], options);
[un_Green,  Emin_Green]  = fmincon(internal_energy_Green_strain, un_natural, [], [], A', b, [], [], [], options);

figure('position', [0, 0, 800, 800]);
set(0,'defaultAxesFontSize',32);
set(0,'defaultTextFontSize',32);

clf;
elastic.draw_individual_rectangles( reshape(un_natural, [2,npoints]) );
xlim([0,12]); ylim([-10,2]); xticks([0:2:12]); grid on;
pbaspect([1 1 1]);
saveas(gcf, 'beam_rectangle_elements_natural_rectangle_elements.png');

clf;
elastic.draw_individual_rectangles( reshape(un_Cauchy, [2,npoints]) );
xlim([0,12]); ylim([-10,2]); xticks([0:2:12]); grid on;
pbaspect([1 1 1]);
saveas(gcf, 'beam_rectangle_elements_deformed_external_force_Cauchy_strain.png');

clf;
elastic.draw_individual_rectangles( reshape(un_Green, [2,npoints]) );
xlim([0,12]); ylim([-10,2]); xticks([0:2:12]); grid on;
pbaspect([1 1 1]);
saveas(gcf, 'beam_rectangle_elements_deformed_external_force_Green_strain.png');

clf;
elastic.draw_individual_rectangles( reshape(un_natural, [2,npoints]) );
elastic.draw_individual_rectangles( reshape(un_Cauchy, [2,npoints]) );
xlim([0,12]); ylim([-10,2]); xticks([0:2:12]); grid on;
pbaspect([1 1 1]);
saveas(gcf, 'beam_rectangle_elements_deformation_external_force_Cauchy.png');

clf;
elastic.draw_individual_rectangles( reshape(un_natural, [2,npoints]) );
elastic.draw_individual_rectangles( reshape(un_Green, [2,npoints]) );
xlim([0,12]); ylim([-10,2]); xticks([0:2:12]); grid on;
pbaspect([1 1 1]);
saveas(gcf, 'beam_rectangle_elements_deformation_external_force_Green.png');

function energy = internal_energy_params( un, npoints, elastic, force )
    disps = reshape(un, [2,npoints]);
    energy = elastic.total_strain_potential_energy(disps);
    energy = energy - force'*disps(:,22);
end

function energy = internal_energy_Green_strain_params( un, npoints, elastic, force )
    disps = reshape(un, [2,npoints]);
    energy = elastic.total_strain_potential_energy_Green_strain(disps);
    energy = energy - force'*disps(:,22);
end
