addpath('../two_dim_fea');

m = 11; n = 3;
width = 10; height = 2; thickness = 1;

[ points, triangles, rectangles ] = rectangular_object( m, n, width, height );
points = points + [ -width/2; 0 ];

% E = 0.1 MPa; rho = 0.020 g/cm^2
Young = 1.0*1e+6; nu = 0.48; density = 0.020;
%pressure = 0.5*1e+7; filename = 'curling_by_inner_pressure/05.png';
%pressure = 1.0*1e+7; filename = 'curling_by_inner_pressure/1.png';
%pressure = 2.0*1e+7; filename = 'curling_by_inner_pressure/2.png';
pressure = 3.0*1e+7; filename = 'curling_by_inner_pressure/3.png';

[ 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);

elastic = elastic.define_subregion([], [1:m-1]);
elastic = elastic.subregion_color( [0.95 0.95 0.95] );
elastic = elastic.define_subregion([], [m:2*(m-1)]);
elastic = elastic.subregion_color( [0.85 0.85 0.85] );

un_natural = zeros(2,npoints);

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

internal_energy_Green_strain = @(un) internal_energy_Green_strain_params( un, npoints, elastic, pressure );

A = zeros(2*npoints, 3); b = zeros(3,1);
k = (m+1)/2;     A(2*k-1, 1) = 1; A(2*k, 2) = 1;
k = (m+1)/2 + m; A(2*k-1, 3) = 1;

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

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

clf;
elastic.draw_individual_rectangles( reshape(un_natural, [2,npoints]) );
xlim([-9,9]); ylim([-2,7]);
xticks([-8:2:8]); yticks([-2:2:7]); grid on;
pbaspect([2 1 1]);
saveas(gcf, 'curling_by_inner_pressure/natural.png');

clf;
elastic.draw_individual_rectangles( reshape(un_Green, [2,npoints]) );
xlim([-9,9]); ylim([-2,7]);
xticks([-8:2:8]); yticks([-2:2:7]); grid on;
pbaspect([2 1 1]);
saveas(gcf, filename);

function energy = internal_energy_Green_strain_params( un, npoints, elastic, pressure )
    persistent index_rects;
    if isempty(index_rects)
        index_rects = elastic.SubRegions(1).Index_Rectangles;
    end

    disps = reshape(un, [2,npoints]);
    energy = elastic.total_strain_potential_energy_Green_strain(disps);
    area_increase = 0;
    for k = index_rects
        rect = elastic.Rectangles(k);
        a = rect.area_increase(disps);
        area_increase = area_increase + a;
    end
    energy = energy  - pressure * elastic.Thickness * area_increase;
end
