function [X,tri] = MovMesh_cube2tet(x,y,z)
%
% usage: [X,tri] = MovMesh_cube2tet(x,y,z)
%
% this function converts a cuboid mesh (defined as a tensor product
% by x and y and z) into a tetrahedral mesh (each subcuboid is divided
% into 6 tetrahedra).
%
% x: points in x-coordinate, of size nx-by-1.
% y: points in y-coordinate, of size ny-by-1.
% z: points in z-coordinate, of size nz-by-1.
% X: (output) coordinates of vertices of the mesh, of size Nv-by-d.
% tri: (output) the connectivity of the mesh, of size N-by-(d+1).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2019 Weizhang Huang (whuang@ku.edu)
imax = max(size(x));
jmax = max(size(y));
kmax = max(size(z));
N = (imax-1)*(jmax-1)*(kmax-1);
[xx, yy, zz] = ndgrid(x,y,z);
X = [xx(:), yy(:), zz(:)];
tri = zeros(6*N,4);
% nodes: 1(0,0,0), 2(0,0,1), 3(1,0,0), 4(1,0,1),
% 5(0,1,0), 6(0,1,1), 7(1,1,0), 8(1,1,1)
% range of vertices: imax*( jmax*(k-1)+(j-1) )+i
I1 = find(X(:,1)