% g, cm, sec

addpath('../two_dim_fea');

% E = 0.1 MPa; c = 0.04 kPa s; rho = 1 g/cm^2
Young = 1.0*1e+6; c = 0.4*1e+3; nu = 0.48; density = 1.00;
[ lambda, mu ] = Lame_constants( Young, nu );
% p = 2e+4 = 0.02e+6 --> 0.002 MPa = 2 KPa
pressure = 2e+4; filename = 'expansion_by_pressure_deformed.png';
%pressure = -2e+4; filename = 'expansion_by_pressure_deformed_negative.png';

mrect = 11; nrect = 2;
[ points, triangles ] = rectangular_object( mrect, nrect, 10, 1 );
thickness = 1;

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

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

clf;
elastic.draw;
xlim([0,10]); ylim([-2,3]); grid on;
%text(5, -1, 'chamber', 'FontSize', 24, 'HorizontalAlignment', 'center');
pbaspect([1 0.5 1]);
saveas(gcf, 'expansion_by_pressure_natural.png');

A = elastic.constraint_matrix([1,11,12,22]);
b = zeros(2*4,1);

index_pressure_area = [mrect:-1:1];

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

internal_energy = @(un) internal_energy_params( un, npoints, elastic, pressure, index_pressure_area );
options = optimset('Display','iter', 'MaxFunEvals',3*10^5);
[un_deformed, Emin] = fmincon(internal_energy, un_natural, [], [], A', b, [], [], [], options);

Emin
clf;
elastic.draw( reshape(un_deformed, [2,npoints]) );
xlim([0,10]); ylim([-2,3]); grid on;
pbaspect([1 0.5 1]);
saveas(gcf, filename);

%disps = reshape(un_deformed, [2,npoints]);
%elastic.surrounded_area_fill([mrect:-1:1], disps(:,mrect:-1:1), [0.2,0.2,0.2]);
%drawnow;

function energy = internal_energy_params( un, npoints, elastic, pressure, index_pressure_area )
    disps = reshape(un, [2,npoints]);
    energy = elastic.total_strain_potential_energy(disps);
    area = elastic.surrounded_area(index_pressure_area, disps(:,index_pressure_area));
    volume = elastic.Thickness * area;
    energy = energy - pressure*volume;
end
