clc; clear;

% parametry geometrie
euler_corner = 20;

% parametry sponge zones
sponge_width = 80;

% parametry velikosti elementu vysledne site
euler_size = 2.0;
sponge_size = 15.0;

% nacteni site, na ktere se vytvori background mesh
soubor_GUI = 0;

if soubor_GUI == 1
    [soubor,cesta,~] = uigetfile('.msh');
    soubor_string = strcat(cesta,soubor);
else
    soubor_string = 'Bajko_Problem1_sponge_geometry.msh';
end
[P,~,E] = MshToMatlab2D(soubor_string);

K = size(E,2); % pocet trojuhelniku
hodnoty = ones(3,K);
% napocitani "profilu" velikosti elementu
hodnoty = euler_size * hodnoty;
% zjemneni ve sponge zones
v_1 = E(1,:); v_2 = E(2,:); v_3 = E(3,:);
x = [P(1,v_1); P(1,v_2); P(1,v_3)];
y = [P(2,v_1); P(2,v_2); P(2,v_3)];
beta_stretch = 1.5;
indexes_L = find((x <= -euler_corner) & (abs(y) <= euler_corner));
hodnoty(indexes_L) = hodnoty(indexes_L) + (sponge_size - euler_size) * (abs(x(indexes_L) + euler_corner) / sponge_width) .^ beta_stretch;
indexes_R = find((x >= euler_corner) & (abs(y) <= euler_corner));
hodnoty(indexes_R) = hodnoty(indexes_R) + (sponge_size - euler_size) * (abs(x(indexes_R) - euler_corner) / sponge_width) .^ beta_stretch;
indexes_B = find((y <= -euler_corner) & (abs(x) <= euler_corner));
hodnoty(indexes_B) = hodnoty(indexes_B) + (sponge_size - euler_size) * (abs(y(indexes_B) + euler_corner) / sponge_width) .^ beta_stretch;
indexes_T = find((y >= euler_corner) & (abs(x) <= euler_corner));
hodnoty(indexes_T) = hodnoty(indexes_T) + (sponge_size - euler_size) * (abs(y(indexes_T) - euler_corner) / sponge_width) .^ beta_stretch;
indexes_LT = find((x <= -euler_corner) & (y >= euler_corner));
vectors = [x(indexes_LT) + euler_corner,y(indexes_LT) - euler_corner];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
hodnoty(indexes_LT) = hodnoty(indexes_LT) + (sponge_size - euler_size) * (distances / sponge_width) .^ beta_stretch;
indexes_LB = find((x <= -euler_corner) & (y <= -euler_corner));
vectors = [x(indexes_LB) + euler_corner,y(indexes_LB) + euler_corner];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
hodnoty(indexes_LB) = hodnoty(indexes_LB) + (sponge_size - euler_size) * (distances / sponge_width) .^ beta_stretch;
indexes_RT = find((x >= euler_corner) & (y >= euler_corner));
vectors = [x(indexes_RT) - euler_corner,y(indexes_RT) - euler_corner];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
hodnoty(indexes_RT) = hodnoty(indexes_RT) + (sponge_size - euler_size) * (distances / sponge_width) .^ beta_stretch;
indexes_RB = find((x >= euler_corner) & (y <= -euler_corner));
vectors = [x(indexes_RB) - euler_corner,y(indexes_RB) + euler_corner];
distances = sqrt(vectors(:,1) .^ 2 + vectors(:,2) .^ 2);
hodnoty(indexes_RB) = hodnoty(indexes_RB) + (sponge_size - euler_size) * (distances / sponge_width) .^ beta_stretch;

filename = 'Bajko_Problem1_sponge_background_mesh.pos';
fileID = fopen(filename,'w');

fprintf(fileID,'View \"background mesh\" {\n');
for i = 1:K
    tri = E(:,i);
    % vrcholy trojuhelniku
    fprintf(fileID,'ST(');
    for j = 1:3
        fprintf(fileID,'%0.8f,%0.8f,0',P(1,tri(j)),P(2,tri(j)));
        if j < 3
            fprintf(fileID,',');
        end
    end
    fprintf(fileID,')');
    % hodnoty ve vrcholech
    fprintf(fileID,'{');
    for j = 1:3
        fprintf(fileID,'%0.8f',hodnoty(j,i));
        if j < 3
            fprintf(fileID,',');
        end
    end
    fprintf(fileID,'};\n');
end
fprintf(fileID,'};');

fclose(fileID);