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

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

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

[x,y] = OptimalNodes2D(N); [r,s] = xytors(x,y);

V = VandermondeMatrix2D(N,r,s); invV = inv(V);
M = invV' * invV; invM = V * V';
[D_r,D_s] = DiffMatrices2D(N,r,s,V,invV);

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

%[r_x,r_y,s_x,s_y,J] = GeometricFactors2D(x,y,D_r,D_s);
[r_x,r_y,s_x,s_y,J] = GeometricFactors2D_my(x,y,D_r,D_s);

fmask1 = find(abs(s + 1) < NODETOL); % cisla uzlu na hrane spojujici 1. a 2. vrchol trojuhelniku
fmask2 = find(abs(r + s) < NODETOL); % cisla uzlu na hrane spojujici 2. a 3. vrchol trojuhelniku
fmask3 = find(abs(r + 1) < NODETOL); % cisla uzlu na hrane spojujici 3. a 1. vrchol trojuhelniku
Fmask = [fmask1; fmask2; fmask3]'; % ve sloupcich cisla ulzu na I, ktere jsou na hranach
Fx = x(Fmask(:),:); Fy = y(Fmask(:),:);

LIFT = Lift2D();

%[n_x,n_y,surfJ] = Normals2D();
[n_x,n_y,surfJ] = Normals2D_my();
Fscale = surfJ ./ J(Fmask(:),:);

[EtoE,EtoF] = Connectivity2D();

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

% Computing weak operators
[V_r,V_s] = GradVandermondeMatrices2D(N,r,s);
D_r_weak = (V * V_r') * M; D_s_weak = (V * V_s') * M;