%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-07-04 16:41:48
 % @LastEditTime: 2021-07-10 19:28:39

 % generate point clouds for each material
%--------------------------

dbstop if error
feature accel on
addpath util util_simulation util_bezier util_micro util_mtimesx util_output

% load paras_simplified
% 
% data_path = ['data\dizhi_rho2_', num2str(globalx),'&', num2str(globaly),'&', num2str(globalz),'.mat'];
% load(data_path);
% 
% [gridx, gridy, gridz] = meshgrid(dx:dx:nelx, dx:dx:nely, dx:dx:nelz);
% gridy = fliplr(gridy);

step = 10;

% for i = 1 : (globaly / step)
%     coorx = gridx((i-1)*step+1 : i*step, :, :);
%     coory = gridy((i-1)*step+1 : i*step, :, :);
%     coorz = gridz((i-1)*step+1 : i*step, :, :);
%     xi    = x((i-1)*10+1 : i*step, :, :);
% 
%     for mat = 1:4
%         idx = find(xi==mat);
%         if ~isempty(idx)
%             cx = coorx(idx);
%             cy = coory(idx);
%             cz = coorz(idx);
% 
%             coor = [cx(:), cy(:), cz(:)];
%             data_path = ['data/point-clouds/', num2str(i),'_', num2str(mat),'.xyz'];
%             fid = fopen(data_path, 'wt') ;
%             fprintf(fid, '%i\n', size(coor, 1)) ;
%             fprintf(fid, '%f %f %f\n', coor.') ;
%             fclose(fid);
%         end
%     end
% end

minx = dx; maxx = nelx;
minz = dx; maxz = nelz;

% for i = [42,49,50] 
for i = [36]
    
    xi = x((i-1)*10+1 : i*step, :, :);

    for mat = 1:2
%         maxy = dx*i*step; 
        miny0 = dx*(i-1)*step;

        idx = find(xi==mat);
        if ~isempty(idx)
            [r,~,~] = ind2sub(size(xi), idx);
            r1 = unique(r);
            
            for j = 1:length(r1)
                r2 = r1(j);

                miny = miny0 + (r2-1) * dx;
                maxy = miny0 + (r2) * dx;

                p = zeros(8,3);
                p(1,:) = [minx, miny, minz];
                p(2,:) = [maxx, miny, minz];
                p(3,:) = [maxx, maxy, minz];
                p(4,:) = [minx, maxy, minz];
                p(5,:) = [minx, miny, maxz];
                p(6,:) = [maxx, miny, maxz];
                p(7,:) = [maxx, maxy, maxz];
                p(8,:) = [minx, maxy, maxz];

                faces = [5, 1, 2, 3, 4;
                        5, 1, 2, 6, 5;
                        5, 2, 3, 7, 6;
                        5, 1, 4, 8, 5;
                        5, 5, 6, 7, 8;
                        5, 3, 4, 8, 7] - 1;

                data_path = ['data/point-clouds/mesh/', num2str(i),'_', num2str(mat),'_part_',num2str(j),'.off'];
                fid = fopen(data_path, 'wt') ;
                fprintf(fid, 'OFF\n');
                fprintf(fid, '8 6 0\n');
                fprintf(fid, '%f %f %f\n', p.') ;
                fprintf(fid, '%i %i %i %i %i\n', faces.') ;
                fclose(fid);
            end
        end
    end
end

for i = [41, 43:48, 51:60]
    maxy = dx*i*step; miny = dx*(i-1)*step;
    
    xi = x((i-1)*10+1 : i*step, :, :);
    for mat = 1:2
        if xi(1)==mat
            p = zeros(8,3);
            p(1,:) = [minx, miny, minz];
            p(2,:) = [maxx, miny, minz];
            p(3,:) = [maxx, maxy, minz];
            p(4,:) = [minx, maxy, minz];
            p(5,:) = [minx, miny, maxz];
            p(6,:) = [maxx, miny, maxz];
            p(7,:) = [maxx, maxy, maxz];
            p(8,:) = [minx, maxy, maxz];

            faces = [5, 1, 2, 3, 4;
                    5, 1, 2, 6, 5;
                    5, 2, 3, 7, 6;
                    5, 1, 4, 8, 5;
                    5, 5, 6, 7, 8;
                    5, 3, 4, 8, 7] - 1;

            data_path = ['data/point-clouds/mesh/', num2str(i),'_', num2str(mat),'.off'];
            fid = fopen(data_path, 'wt') ;
            fprintf(fid, 'OFF\n');
            fprintf(fid, '8 6 0\n');
            fprintf(fid, '%f %f %f\n', p.') ;
            fprintf(fid, '%i %i %i %i %i\n', faces.') ;
            fclose(fid);
        end
    end
end

% %% material 3 & 4
% for k = 1:10
    
%     flags5 = false(globalx, globalz);
%     flags6 = false(globalx, globalz);
%     for i = 1:globalx
%         for j = 1:globalz                
            
%             if (i - rand_x3(k))^2/rand_a3(k)^2 + (j - rand_z3(k))^2/rand_b3(k)^2 < 1
%                 flags5(i, j) = 1;
%             end                        
%             if (i - rand_x4(k))^2/rand_a4(k)^2 + (j - rand_z4(k))^2/rand_b4(k)^2 < 1
%                 flags6(i, j) = 1;
%             end            
%         end
%     end
%     cx3 = gridx(rand_y3(k):rand_y3(k)+thichkness_y, flags5);
%     cy3 = gridy(rand_y3(k):rand_y3(k)+thichkness_y, flags5);
%     cz3 = gridz(rand_y3(k):rand_y3(k)+thichkness_y, flags5);

%     cx4 = gridx(rand_y4(k):rand_y4(k)+thichkness_y, flags6);
%     cy4 = gridy(rand_y4(k):rand_y4(k)+thichkness_y, flags6);
%     cz4 = gridz(rand_y4(k):rand_y4(k)+thichkness_y, flags6);
    
%     point_clouds3 = zeros(length(cx3(:)), 3);
%     parfor ii = 1:length(cx3(:))
%         cx3_1 = cx3(ii);
%         cy3_1 = cy3(ii);
%         cz3_1 = cz3(ii);

%         point_clouds3(ii, :) = [cx3_1, cy3_1, cz3_1];
%     end
%     data_path = ['data/stl/', num2str(k),'_', num2str(3),'.stl'];
%     fid = fopen(data_path, 'wt') ;
%     fprintf(fid, '%f %f %f\n', point_clouds3.') ;
%     fclose(fid);

%     point_clouds4 = zeros(length(cx4(:)), 3);
%     for ii = 1:length(cx4(:))
%         cx4_1 = cx4(ii);
%         cy4_1 = cy4(ii);
%         cz4_1 = cz4(ii);

%         point_clouds4(ii, :) = [cx4_1, cy4_1, cz4_1];
%     end  
%     data_path = ['data/stl/', num2str(k),'_', num2str(4),'.stl'];
%     fid = fopen(data_path, 'wt') ;
%     fprintf(fid, '%f %f %f\n', point_clouds4.') ;
%     fclose(fid);
% end

% %% generate stl
% nelx = 42;  
% nely = 20;
% nelz = 30;
% microx = 30;
% globaly = nely * microx;
% globalx = nelx * microx;
% globalz = nelz * microx;

% [gridx, gridy, gridz] = meshgrid(1:globalx, 1:globaly, 1:globalz);
% gridy = fliplr(gridy);

% load('paras_modeling.mat', 'x');
% for mat = 1:4
%     fv = isosurface(gridx, gridy, gridz, x, mat, 'verbose');
%     if ~isempty(fv.vertices)
%         data_path = ['data/stl/mat', num2str(mat),'.stl'];
%         stlwrite(data_path, fv);
%     end
% end

% for mat = 1:4
%     for i = 1:60
%         cx = gridx((i-1)*10+1:i*10,:,:);
%         cy = gridy((i-1)*10+1:i*10,:,:);
%         cz = gridz((i-1)*10+1:i*10,:,:);
%         xi = x((i-1)*10+1:i*10,:,:);
        
%         fv = isosurface(cx, cy, cz, xi, mat, 'verbose');
%         if ~isempty(fv.vertices)
%             data_path = ['data/stl/', num2str(i),'_', num2str(mat),'.stl'];
%             stlwrite(data_path, fv);
%         end
%     end
% end