clc;
clear;
close all;

N_faces = 4;

mesh_plot = 1;
pause_plot = 0;

N = 1;
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);

method_string = 'Gmsh - Delaunay and Optimize';
[P,E] = MshToMatlab3D('Acoustic_pulse_propagation_small_h=2.00_Delaunay.msh');
MeshQualityAnalysis;
if mesh_plot == 1
    figure;
    tetramesh(E',P');
    axis equal;
    xlim([-5, 15]); ylim([-5, 15]); zlim([-5, 15]);
end

method_string = 'Gmsh - Delaunay and Optimize (Netgen)';
[P,E] = MshToMatlab3D('Acoustic_pulse_propagation_small_h=2.00_Delaunay_Netgen.msh');
MeshQualityAnalysis;
if mesh_plot == 1
    figure;
    tetramesh(E',P');
    axis equal;
    xlim([-5, 15]); ylim([-5, 15]); zlim([-5, 15]);
end

%method_string = 'FLUENT Gambit';
%[Nv, VX, VY, VZ, K, EToV] = MeshReaderGambit3D('advection_3D_100.neu');
%P = [VX; VY; VZ];
%E = EToV';
%MeshQualityAnalysis;
%if mesh_plot == 1
%    figure;
%    tetramesh(E',P');
%    axis equal;
%    xlim([-5, 15]); ylim([-5, 15]); zlim([-5, 15]);
%end