clc;

Globals3D;

N_faces = 4;
NODETOL = 1e-9;

N_vect = 1:6;
n = length(N_vect);

format_spec = '%0.15f';
for i = 1:n
    N = N_vect(i);
    N_p = ((N + 1) * (N + 2) * (N + 3)) / 6;
    N_fp = ((N + 1) * (N + 2)) / 2;
    
    % zatim jsou r, s radkove vektory, proto ukladam transpozici
    x_opt_eq = load(sprintf('x_opt_eq_Hesthaven_N=%d.txt',N));
    y_opt_eq = load(sprintf('y_opt_eq_Hesthaven_N=%d.txt',N));
    z_opt_eq = load(sprintf('z_opt_eq_Hesthaven_N=%d.txt',N));
    [r,s,t] = xyztorst(x_opt_eq,y_opt_eq,z_opt_eq);
    filename = sprintf('r_3D_N=%d.txt',N);
    MatrixToTextFile(r',filename,format_spec,' ');
    filename = sprintf('s_3D_N=%d.txt',N);
    MatrixToTextFile(s',filename,format_spec,' ');
    filename = sprintf('t_3D_N=%d.txt',N);
    MatrixToTextFile(t',filename,format_spec,' ');
    
    V = VandermondeMatrix3D(N,r,s,t); invV = inv(V);
    filename = sprintf('V_3D_N=%d.txt',N);
    MatrixToTextFile(V,filename,format_spec,' ');
    filename = sprintf('invV_3D_N=%d.txt',N);
    MatrixToTextFile(invV,filename,format_spec,' ');
    
    M = invV' * invV; invM = V * V';
    filename = sprintf('M_3D_N=%d.txt',N);
    MatrixToTextFile(M,filename,format_spec,' ');
    filename = sprintf('invM_3D_N=%d.txt',N);
    MatrixToTextFile(invM,filename,format_spec,' ');
    
    V_r = load(sprintf('V_r_Hesthaven_N=%d.txt',N));
    V_s = load(sprintf('V_s_Hesthaven_N=%d.txt',N));
    V_t = load(sprintf('V_t_Hesthaven_N=%d.txt',N));
    filename = sprintf('V_r_3D_N=%d.txt',N);
    MatrixToTextFile(V_r,filename,format_spec,' ');
    filename = sprintf('V_s_3D_N=%d.txt',N);
    MatrixToTextFile(V_s,filename,format_spec,' ');
    filename = sprintf('V_t_3D_N=%d.txt',N);
    MatrixToTextFile(V_t,filename,format_spec,' ');
    
    [D_r,D_s,D_t] = DiffMatrices3D(N,r,s,t,V,invV);
    filename = sprintf('D_r_3D_N=%d.txt',N);
    MatrixToTextFile(D_r,filename,format_spec,' ');
    filename = sprintf('D_s_3D_N=%d.txt',N);
    MatrixToTextFile(D_s,filename,format_spec,' ');
    filename = sprintf('D_t_3D_N=%d.txt',N);
    MatrixToTextFile(D_t,filename,format_spec,' ');
    
    D_r_weak = (V * V_r') * M; D_s_weak = (V * V_s') * M; D_t_weak = (V * V_t') * M;
    filename = sprintf('D_r_weak_3D_N=%d.txt',N);
    MatrixToTextFile(D_r_weak,filename,format_spec,' ');
    filename = sprintf('D_s_weak_3D_N=%d.txt',N);
    MatrixToTextFile(D_s_weak,filename,format_spec,' ');
    filename = sprintf('D_t_weak_3D_N=%d.txt',N);
    MatrixToTextFile(D_t_weak,filename,format_spec,' ');
    
    fmask1 = find(abs(t + 1) < NODETOL);
    fmask2 = find(abs(s + 1) < NODETOL);
    fmask3 = find(abs(r + s + t + 1) < NODETOL);
    fmask4 = find(abs(r + 1) < NODETOL);
    Fmask = [fmask1; fmask2; fmask3; fmask4]';
    filename = sprintf('Fmask_3D_N=%d.txt',N);
    MatrixToTextFile(Fmask,filename,'%d',' ');

    LIFT = Lift3D();
    filename = sprintf('LIFT_3D_N=%d.txt',N);
    MatrixToTextFile(LIFT,filename,format_spec,' ');
end