clc; clear;

floatformat = '%0.2f';

Globals2D;

fileID = fopen('AlgorithmDemonstration2D.txt','w');

NODETOL = 1e-12;
N_faces = 3;


disp('Order of local approximations:');
N = 3
disp('');
N_fp = N + 1;
N_p = ((N + 1) * (N + 2)) / 2;

disp('Number of points within each element:');
N_p
disp('');
disp('Optimized points on reference element:');
[x,y] = OptimalNodes2D(N);
[r,s] = xytors(x,y);
disp(x);
disp('');
disp(y);
disp('');
disp(r);
disp('');
disp(s);
disp('');
fprintf(fileID,'x:\n');
LatexCodeMatrix(x,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'y:\n');
LatexCodeMatrix(y,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'r:\n');
LatexCodeMatrix(r,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'s:\n');
LatexCodeMatrix(s,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Vandermonde matrix:');
V = VandermondeMatrix2D(N,r,s); invV = inv(V);
disp(V);
disp('');
fprintf(fileID,'V:\n');
LatexCodeMatrix(V,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Mass matrix M on the reference element:');
M = inv(V * V'); invM = V * V';
disp(M);
disp('');
fprintf(fileID,'M:\n');
LatexCodeMatrix(M,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Vandermonde matrices of derivatives V_r and V_s:');
[V_r,V_s] = GradVandermondeMatrices2D(N,r,s);
disp(V_r);
disp('');
disp(V_s);
disp('');
fprintf(fileID,'V_r:\n');
LatexCodeMatrix(V_r,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'V_s:\n');
LatexCodeMatrix(V_s,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Differentiation matrices D_r and D_s:');
[D_r,D_s] = DiffMatrices2D(N,r,s,V,invV);
disp(D_r);
disp('');
disp(D_s);
disp('');
fprintf(fileID,'D_r:\n');
LatexCodeMatrix(D_r,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'D_s:\n');
LatexCodeMatrix(D_s,floatformat,fileID);
fprintf(fileID,'\n\n');


disp('Mesh for demonstration:');
[P,B,E] = MshToMatlab2D('Advection2D_4elem.msh');
disp(P);
disp('');
disp(E);
disp('');
disp('Number of elements:');
K = size(E,2)
fprintf(fileID,'P:\n');
LatexCodeMatrix(P,'%0.1f',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'E:\n');
LatexCodeMatrix(E,'%d',fileID);
fprintf(fileID,'\n\n');

disp('Coordinates of all points in the grid:');
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));
disp(x);
disp('');
disp(y);
disp('');
fprintf(fileID,'x:\n');
LatexCodeMatrix(x,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'y:\n');
LatexCodeMatrix(y,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Geometric factors; Jacobians:');
[r_x,r_y,s_x,s_y,J] = GeometricFactors2D(x,y,D_r,D_s);
disp(J);
disp('');
disp(r_x)
disp('');
disp(r_y)
disp('');
disp(s_x)
disp('');
disp(s_y)
disp('');
fprintf(fileID,'J:\n');
LatexCodeMatrix(J,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'r_x:\n');
LatexCodeMatrix(r_x,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'r_y:\n');
LatexCodeMatrix(r_y,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'s_x:\n');
LatexCodeMatrix(s_x,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'s_y:\n');
LatexCodeMatrix(s_y,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Extracting information on faces; local numbers of nodes on faces');
fmask1 = find(abs(s + 1) < NODETOL);
fmask2 = find(abs(r + s) < NODETOL);
fmask3 = find(abs(r + 1) < NODETOL);
Fmask = [fmask1; fmask2; fmask3]';
disp(Fmask);
disp('');
disp('Coordinates of all points on faces:');
Fx = x(Fmask(:),:); Fy = y(Fmask(:),:);
disp(Fx);
disp('');
disp(Fy);
disp('');
fprintf(fileID,'Fmask:\n');
LatexCodeMatrix(Fmask,'%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'Fx:\n');
LatexCodeMatrix(Fx,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'Fy:\n');
LatexCodeMatrix(Fy,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Operator for integrals along boundary:');
LIFT = Lift2D();
disp(LIFT);
disp('');
fprintf(fileID,'LIFT:\n');
LatexCodeMatrix(LIFT,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Normals in points on faces:');
[n_x,n_y,surfJ] = Normals2D();
disp(n_x);
disp('');
disp(n_y);
disp('');
fprintf(fileID,'n_x:\n');
LatexCodeMatrix(n_x,floatformat,fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'n_y:\n');
LatexCodeMatrix(n_y,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Inverse Jacobians in points on faces:');
Fscale = surfJ ./ J(Fmask(:),:);
disp(Fscale);
disp('');
fprintf(fileID,'Fscale:\n');
LatexCodeMatrix(Fscale,floatformat,fileID);
fprintf(fileID,'\n\n');

disp('Connectivity of the grid; element to element connections:');
[EtoE,EtoF] = Connectivity2D(E);
disp(full(EtoE));
disp('');
disp('element to face connections:');
disp(full(EtoF));
disp('');
fprintf(fileID,'FtoV:\n');
LatexCodeMatrix(full(spFtoV),'%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'FtoV * FtoV^T:\n');
LatexCodeMatrix(full(spFtoV * spFtoV'),'%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'FtoV * FtoV^T - 2I:\n');
LatexCodeMatrix(full(spFtoF),'%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'EtoE:\n');
LatexCodeMatrix(EtoE,'%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'EtoF:\n');
LatexCodeMatrix(EtoF,'%d',fileID);
fprintf(fileID,'\n\n');

disp('Vertex connectivity and boundary; connectivity of vertices:');
[mapM, vmapM, mapP, vmapP, mapB, vmapB] = BuildMaps2D();
[mapM'; mapP']
[vmapM'; vmapP']
disp('Local numbering of boundary nodes:');
disp(mapB');
disp('');
disp('Global numbering of boundary nodes:');
disp(vmapB');
disp('');
fprintf(fileID,'mapM:\n');
LatexCodeMatrix(mapM','%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'vmapM:\n');
LatexCodeMatrix(vmapM','%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'mapP:\n');
LatexCodeMatrix(mapP','%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'vmapP:\n');
LatexCodeMatrix(vmapP','%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'mapB:\n');
LatexCodeMatrix(mapB','%d',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'vmapB:\n');
LatexCodeMatrix(vmapB','%d',fileID);

fclose(fileID);