%This script will plot a given mesh file in .ncdf format generated from the
%acdtool meshconvert routine. The mesh file can be of type 'full', 'halfx',
%'halfy', or 'quarter' if the mesh contains planes of symmetry. In such
%cases the mesh is tiled via reflection across the yz-plane, xz-plane, or
%both for display purposes.
%
%Written by David Bizzozero
%March 31, 2021

% mesh_file = 'quarter_mesh_demo.ncdf';
% type = 'quarter';

mesh_file = 'SRF-GUN-cathR10mm.ncdf';  %Filename for .ncdf mesh
type = 'quarter';                      %Mesh symmetry type

%% End of user inputs section

%Load interior tetrahedra
t_ext = ncread(mesh_file,'tetrahedron_exterior');
C = ncread(mesh_file,'coords');
X = C(1,:)';    Y = C(2,:)';    Z = C(3,:)';

Tets = t_ext(2:5,:)'+1;
Tris = zeros(4*size(Tets,1),3);

%Create 4 triangles from each boundary tetrahedron
for i = 1:size(Tets,1)
    Tris(4*i-3,:) = [Tets(i,1) Tets(i,3) Tets(i,2)];
    Tris(4*i-2,:) = [Tets(i,1) Tets(i,4) Tets(i,3)];
    Tris(4*i-1,:) = [Tets(i,1) Tets(i,2) Tets(i,4)];
    Tris(4*i,:)   = [Tets(i,2) Tets(i,3) Tets(i,4)];
end

%Locate boundary triangles on from boundary tetrahedra (boundary code 6)
Tflag = t_ext(6:9,:)==6;
Tris = Tris(Tflag(:),:);

%Mirror boundary triangles using symmetry if needed
if strcmpi(type,'halfx')
    verts = size(X,1);
    X = [X;X];
    Y = [Y;-Y];
    Z = [Z;Z];
    Tris = [Tris;Tris+verts];
elseif strcmpi(type,'halfy')
    verts = size(X,1);
    X = [X;-X];
    Y = [Y;Y];
    Z = [Z;Z];
    Tris = [Tris;Tris+verts];
elseif strcmpi(type,'quarter')
    verts = size(X,1);
    X = [X;X;-X;-X];
    Y = [Y;-Y;-Y;Y];
    Z = [Z;Z;Z;Z];
    Tris = [Tris;Tris+verts;Tris+2*verts;Tris+3*verts];
end

%Plot tetrahedra and Cartesian grid near z=0 and view from front
fontsz = 20;
figure('position',get(0,'screensize'))
tm = trimesh(Tris,Z,X,Y,ones(size(X)),...
    'facealpha',1,'edgecolor','k','facecolor','c');
hold on
axis equal
set(gca,'fontsize',fontsz,'fontname','times')
xlabel('$z\,{\rm [m]}$','interpreter','latex','fontsize',fontsz)
ylabel('$x\,{\rm [m]}$','interpreter','latex','fontsize',fontsz)
zlabel('$y\,{\rm [m]}$','interpreter','latex','fontsize',fontsz)