rk4_ls;

NODETOL = 1e-12;

N_v = size(P,2);
K = size(E,2);

%[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_my(x,y,D_r,D_s);
%
%fmask1 = find(abs(s + 1) < NODETOL);
%fmask2 = find(abs(r + s) < NODETOL);
%fmask3 = find(abs(r + 1) < NODETOL);
%Fmask = [fmask1; fmask2; fmask3]'; 
%Fx = x(Fmask(:),:); Fy = y(Fmask(:),:);
%
%[n_x,n_y,surfJ] = Normals2D_my();
%Fscale = surfJ ./ J(Fmask(:),:);

[EtoE,EtoF] = Connectivity2D();

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