clc; clear;

Globals1D;

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

NODETOL = 1e-12;
N_faces = 2;
N_fp = 1;

disp('Order of local approximations:');
N = 4
disp('');
N_p = N + 1;

disp('Generated normalized Legendre polynomials:');
LegPols = LegendrePolynomialsNormalized(N);
modal_basis = zeros(N_p);
for i = 1:N_p
    modal_basis(i,(N_p - (i - 1)):N_p) = LegPols{i};
end
disp(modal_basis);
disp('');
fprintf(fileID,'modal basis:\n');
for i = 1:N_p
    LatexCodePolynomial(modal_basis(i,:),'r',fileID);
    fprintf(fileID,'\n');
end
fprintf(fileID,'\n');

disp('The simple basis:');
polynomials = zeros(N + 1);
for i = 1:(N + 1)
    polynomials(i,(N + 1) - i + 1) = 1;
end
disp(polynomials);
disp('');
fprintf(fileID,'simple basis:\n');
for i = 1:N_p
    LatexCodePolynomial(polynomials(i,:),'r',fileID);
    fprintf(fileID,'\n');
end
fprintf(fileID,'\n');

disp('Simple basis after Gram-Schmidt orthonormalization:');
orthonormalized = PolynomialGramSchmidtOrthonormalization(polynomials);
disp(orthonormalized);
disp('');
fprintf(fileID,'orthonormalized simple basis:\n');
for i = 1:N_p
    LatexCodePolynomial(orthonormalized(i,:),'r',fileID);
    fprintf(fileID,'\n');
end
fprintf(fileID,'\n');

disp('Number of points within each element:');
N_p
disp('');
disp('Legendre-Gauss-Lobatto nodes on reference element:');
r = GaussLobattoNodes(N_p);
disp(r);
disp('');
fprintf(fileID,'r:\n');
LatexCodeMatrix(r,'%f',fileID);
fprintf(fileID,'\n\n');

disp('Vandermonde matrix:');
V = VandermondeMatrix1D(N,r);
disp(V);
disp('');
fprintf(fileID,'V:\n');
LatexCodeMatrix(V,'%f',fileID);
fprintf(fileID,'\n\n');

disp('Nodal basis after transformation from modal basis:');
nodal_basis = inv(V') * modal_basis;
disp(nodal_basis);
disp('');
fprintf(fileID,'nodal basis:\n');
for i = 1:N_p
    LatexCodePolynomial(nodal_basis(i,:),'r',fileID);
    fprintf(fileID,'\n');
end
fprintf(fileID,'\n');

disp('Lagrange polynomials defined by nodes r:');
LagPols = LagrangePolynomialBasis(r);
LagPols_matrix = zeros(N_p);
for i = 1:N_p
    LagPols_matrix(i,:) = LagPols{i};
end
disp(LagPols_matrix);
disp('');
fprintf(fileID,'Lagrange polynomials:\n');
for i = 1:N_p
    LatexCodePolynomial(LagPols_matrix(i,:),'r',fileID);
    fprintf(fileID,'\n');
end
fprintf(fileID,'\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,'%f',fileID);
fprintf(fileID,'\n\n');

disp('Differentiation matrix D_r:');
D_r = DiffMatrix1D(N,r,V);
disp(D_r);
disp('');
fprintf(fileID,'D_r:\n');
LatexCodeMatrix(D_r,'%f',fileID);
fprintf(fileID,'\n\n');

disp('Operator for integrals along boundary');
LIFT = Lift1D();
disp(LIFT);
disp('');
fprintf(fileID,'LIFT:\n');
LatexCodeMatrix(LIFT,'%f',fileID);
fprintf(fileID,'\n\n');


disp('Mesh for demonstration:');
P = [0.0, 1.0, 2.0, 3.0]
E = [1, 2, 3; 2, 3, 4]
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,:);
x = ones(N_p,1) * P(v_1) + 0.5 * (r + 1)' * (P(v_2) - P(v_1));
disp(x);
disp('');
fprintf(fileID,'x:\n');
LatexCodeMatrix(x,'%f',fileID);
fprintf(fileID,'\n\n');

disp('Geometric factors; Jacobians:');
[r_x,J] = GeometricFactors1D(x,D_r);
disp(J);
disp('');
disp('and inverse Jacobians:');
disp(r_x)
disp('');
fprintf(fileID,'J:\n');
LatexCodeMatrix(J,'%f',fileID);
fprintf(fileID,'\n\n');
fprintf(fileID,'r_x:\n');
LatexCodeMatrix(r_x,'%f',fileID);
fprintf(fileID,'\n\n');

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

disp('Normals in points on faces:');
n_x = Normals1D();
disp(n_x);
disp('');
fprintf(fileID,'n_x:\n');
LatexCodeMatrix(n_x,'%1.1f',fileID);
fprintf(fileID,'\n\n');

disp('Inverse Jacobians in points on faces:');
Fscale = r_x(Fmask,:);
disp(Fscale);
disp('');
fprintf(fileID,'Fscale:\n');
LatexCodeMatrix(Fscale,'%f',fileID);
fprintf(fileID,'\n\n');

disp('Connectivity of the grid; element to element connections:');
[EtoE,EtoF] = Connectivity1D(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 - I:\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:');
[vmapM,vmapP,mapB,vmapB] = BuildMaps1D();
[vmapM'; vmapP']
disp('Local numbering of boundary nodes:');
disp(mapB');
disp('');
disp('Global numbering of boundary nodes:');
disp(vmapB');
disp('');
fprintf(fileID,'vmapM:\n');
LatexCodeMatrix(vmapM','%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);