% Comparing calculations based on Cauchy and Green strains
% left end fixed and gravitational force acts downward

% コーシー歪みとグリーン歪みによる計算を比較
% ビームの左端を固定．下向きに重力が作用する．

% g, cm, sec

addpath('../two_dim_fea');

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

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

% E = 0.1 MPa; rho = 10.00 g/cm^2
Young = 1.0*1e+6; nu = 0.48; density = 10.00;
[ lambda, mu ] = Lame_constants( Young, nu );
% g = 9.8 m/s^2 = 980 cm/s^2
grav = [0; -980];

npoints = size(points,2);
ntriangles = size(triangles,1);
elastic = Body(npoints, points, ntriangles, triangles, thickness);
elastic = elastic.mechanical_parameters(density, lambda, mu);

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

strain_and_gravitational_energy = @(un) strain_and_gravitational_energy_params( un, npoints, elastic, grav );
strain_and_gravitational_energy_Green_strain = @(un) strain_and_gravitational_energy_Green_strain_params( un, npoints, elastic, grav );

% 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_g, Emin_Cauchy_g] = fmincon(strain_and_gravitational_energy, un_natural, [], [], A', b, [], [], [], options);
[un_Green_g,  Emin_Green_g]  = fmincon(strain_and_gravitational_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( reshape(un_natural, [2,npoints]) );
xlim([0,12]); ylim([-10,2]); xticks([0:2:12]); grid on;
pbaspect([1 1 1]);
saveas(gcf, 'beam_natural.png');

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

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

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

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

function energy = strain_and_gravitational_energy_params( un, npoints, elastic, grav )
    disps = reshape(un, [2,npoints]);
    energy = elastic.total_strain_potential_energy(disps) ...
           + elastic.total_gravitational_potential_energy(disps, grav);
end

function energy = strain_and_gravitational_energy_Green_strain_params( un, npoints, elastic, grav )
    disps = reshape(un, [2,npoints]);
    energy = elastic.total_strain_potential_energy_Green_strain(disps) ...
           + elastic.total_gravitational_potential_energy(disps, grav);
end
