topo_file_name = 'VIDAL_CA-24000.dem';
[latgrat, longgrat, topo_mat, dem_info]   = usgs24kdem(topo_file_name, 1);
[topo_mat_rows, topo_mat_cols]            = size(topo_mat)
x_vec                                     = 30 * linspace(0, topo_mat_cols - 1, topo_mat_cols);
y_vec                                     = -30 * linspace(0, topo_mat_rows - 1, topo_mat_rows);
[x_mat, y_mat]                            = meshgrid(x_vec, y_vec);
node = zeros(8*length(y_vec)*length(x_vec),4);
face = zeros(6*length(y_vec)*length(x_vec),5);
for y=1:length(y_vec)
    for x=1:length(x_vec)
        blocknum=length(x_vec)*(y-1)+x-1;
        node(8*blocknum+1,:) = [8*blocknum+1,x_vec(x)-15,y_vec(y)-15,topo_mat(y,x)];
        node(8*blocknum+2,:) = [8*blocknum+2,x_vec(x)+15,y_vec(y)-15,topo_mat(y,x)];
        node(8*blocknum+3,:) = [8*blocknum+3,x_vec(x)+15,y_vec(y)+15,topo_mat(y,x)];
        node(8*blocknum+4,:) = [8*blocknum+4,x_vec(x)-15,y_vec(y)+15,topo_mat(y,x)];
        node(8*blocknum+5,:) = [8*blocknum+5,x_vec(x)-15,y_vec(y)-15,0];
        node(8*blocknum+6,:) = [8*blocknum+6,x_vec(x)+15,y_vec(y)-15,0];
        node(8*blocknum+7,:) = [8*blocknum+7,x_vec(x)+15,y_vec(y)+15,0];
        node(8*blocknum+8,:) = [8*blocknum+8,x_vec(x)-15,y_vec(y)+15,0];
        face(6*blocknum+1,:) = [1,1,2,3,4]+6*blocknum;
        face(6*blocknum+2,:) = [2,8,7,6,5]+6*blocknum;
        face(6*blocknum+3,:) = [3,8,5,1,4]+6*blocknum;
        face(6*blocknum+4,:) = [4,6,5,3,2]+6*blocknum;
        face(6*blocknum+5,:) = [5,5,6,2,1]+6*blocknum;
        face(6*blocknum+6,:) = [6,7,8,4,3]+6*blocknum;
    end
end
[x_mat, y_mat] = meshgrid(x_vec, y_vec);