rk4_ls; % definition of coefficients for low-storage variant of Runge-Kutta method

% lokalni reprezentace reseni
N_p = ((N + 1) * (N + 2) * (N + 3)) / 6; % pocet uzlu na 1 elementu = pocet neznamych parametru (koeficientu interpol. polynomu na elementu)
N_fp = ((N + 1) * (N + 2)) / 2; % pocet uzlu na 1 stene
N_faces = 4; % 4 steny pro kazdy element ve 3D
NODETOL = 1e-9;

N_v = size(P,2); % pocet uzlu v siti
K = size(E,2); % pocet elementu v siti

%[x_opt_eq,y_opt_eq,z_opt_eq] = OptimalNodes3D(N); [r,s,t] = xyztorst(x,y,z);
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);

V = VandermondeMatrix3D(N,r,s,t); invV = inv(V);
M = invV' * invV; invM = V * V';
[D_r,D_s,D_t] = DiffMatrices3D(N,r,s,t,V,invV);

v_1 = E(1,:); v_2 = E(2,:); v_3 = E(3,:); v_4 = E(4,:);
x = 0.5 * (-(r + s + t + 1)' * P(1,v_1) + (r + 1)' * P(1,v_2) + (s + 1)' * P(1,v_3) + (t + 1)' * P(1,v_4));
y = 0.5 * (-(r + s + t + 1)' * P(2,v_1) + (r + 1)' * P(2,v_2) + (s + 1)' * P(2,v_3) + (t + 1)' * P(2,v_4));
z = 0.5 * (-(r + s + t + 1)' * P(3,v_1) + (r + 1)' * P(3,v_2) + (s + 1)' * P(3,v_3) + (t + 1)' * P(3,v_4));

[r_x,r_y,r_z,s_x,s_y,s_z,t_x,t_y,t_z,J] = GeometricFactors3D_my(x,y,z,D_r,D_s,D_t);

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]';
Fx = x(Fmask(:),:); Fy = y(Fmask(:),:); Fz = z(Fmask(:),:);

LIFT = Lift3D();

[n_x,n_y,n_z,surfJ] = Normals3D_my();
Fscale = surfJ ./ J(Fmask(:),:);

[EtoE,EtoF] = Connectivity3D();

[mapM, vmapM, mapP, vmapP, mapB, vmapB] = BuildMaps3D();

% Computing weak operators
%[V_r,V_s,V_t] = GradVandermondeMatrices3D(N,r,s,t);
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));
D_r_weak = (V * V_r') * M; D_s_weak = (V * V_s') * M; D_t_weak = (V * V_t') * M;