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 = []
face = []
for y=1:length(y_vec)
    for x=1:length(x_vec)
        blocknum=length(x_vec)*(y-1)+x;
        node = [node;[x_vec(x)-15,y_vec(y)-15,topo_mat(x,y)]];
        node = [node;[x_vec(x)+15,y_vec(y)-15,topo_mat(x,y)]];
        node = [node;[x_vec(x)+15,y_vec(y)+15,topo_mat(x,y)]];
        node = [node;[x_vec(x)-15,y_vec(y)+15,topo_mat(x,y)]];
        node = [node;[x_vec(x)-15,y_vec(y)-15,0]];
        node = [node;[x_vec(x)+15,y_vec(y)-15,0]];
        node = [node;[x_vec(x)+15,y_vec(y)+15,0]];
        node = [node;[x_vec(x)-15,y_vec(y)+15,0]];
        face = [face;[0,1,2,3]+6*blocknum];
        face = [face;[7,6,5,4]+6*blocknum];
        face = [face;[7,4,0,3]+6*blocknum];
        face = [face;[5,6,2,1]+6*blocknum];
        face = [face;[4,5,1,0]+6*blocknum];
        face = [face;[6,7,3,2]+6*blocknum];
    end
end
[x_mat, y_mat] = meshgrid(x_vec, y_vec);
node
face