%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-07-12 09:23:41
 % @LastEditTime: 2021-08-15 21:21:49
%--------------------------
dbstop if error
feature accel on
addpath util util_simulation util_bezier util_micro util_mtimesx util_output

%% adjustable parameters
K0 = 7e6;    G0 = 3.2e6;
K1 = 1.86e7; G1 = 9e6;
K2 = 1.38e8; G2 = 5.96e7;
K3 = 6.3e8;  G3 = 3.86e8;

nelx = 42;  
nely = 20;
nelz = 30;
microx = 30;
lx = 1;
num_add = 0; % additional macro dofs on the border

optDesign = 'bending2'; vF = 1e5; offset = 2;

opt_load_edofmat = 0;
opt_load_rho = 0;

%% init
dx = lx / microx;

[~, ~, edofMat_mi] = prepare_assemble(microx, microx, microx);
[boundarys, inner] = micro_boundary(microx);
newdofs = [boundarys(:); inner(:)];
dofid   = reorder_dofs(newdofs); 

if opt_load_rho==1
    globaly = nely * microx;
    globalx = nelx * microx;
    globalz = nelz * microx;

    data_path = ['data\dizhi_rho2_', num2str(globalx),'&', num2str(globaly),'&', num2str(globalz),'.mat'];
    load(data_path);
end

%% reconstruction
% data_path = ['data/U2_',optDesign,'_F',num2str(vF),'_offset',num2str(offset),'.mat'];
% load(data_path);
% 
% if opt_load_edofmat==1
%     load('dizhi_edofMat.mat', 'edofMat_ma');
% end
% 
% globaly = 3 * microx;
% globalx = nelx * microx;
% globalz = nelz * microx;
% uglobal = zeros(globalx*globaly*globalz, 3);
% 
% [gridx, gridy, gridz] = meshgrid(dx:dx:nelx, dx:dx:3, dx:dx:nelz);
% gridy = fliplr(gridy);

% for i = 1:nelx
%     for j = 1:3
%         data_path2 = ['output_data/sK_new2/R/', num2str(nelx),'&', num2str(nely),'&', num2str(nelz),'_', ...
%             num2str(microx),'m_',num2str(num_add),'pnts_',num2str(i),'_',num2str(j),'.mat'];        
%         load(data_path2); % load regulated_matrix
% 
%         for k = 1:nelz
%             ele  = (k-1)*nelx*nely + (i-1)*nely + j;
%             edof = edofMat_ma(ele, :);
%             ue = U(edof);
%             he = regulated_tmp{k} * ue; % he order: [4-borders; inner]
%             u0 = he(dofid);
%             u1 = u0(edofMat_mi); % [microx^3, 24];
% 
%             umicro(:,1) = mean(u1(:,1:3:end), 2);
%             umicro(:,2) = mean(u1(:,2:3:end), 2);
%             umicro(:,3) = mean(u1(:,3:3:end), 2);            
% 
%             for ii = 1:microx
%                 for jj = 1:microx
%                     for kk = 1:microx
%                         % which micro element in full domain
%                         fulli = (i-1)*microx + ii;
%                         fullj = (j-1)*microx + jj;
%                         fullk = (k-1)*microx + kk;
%                         global_ele = (fullk-1)*globalx*globaly + (fulli-1)*globaly + fullj;
% 
%                         % which micro element in this macro-ele
%                         local_ele = (kk-1)*microx^2 + (ii-1)*microx + jj;
% 
%                         uglobal(global_ele, :) = umicro(local_ele, :);
%                     end
%                 end
%             end
%         end
%     end
% end
% 
% defx = reshape(uglobal(:,1), globaly, globalx, globalz) + gridx;
% defy = reshape(uglobal(:,2), globaly, globalx, globalz) + gridy;
% defz = reshape(uglobal(:,3), globaly, globalx, globalz) + gridz;
% clear gridx gridy gridz uglobal

% load def
% defy0 = defy - gridy;
% gridy0 = flipud(gridy);
% defy = defy0 + gridy0;

% load def_correct_y
cube = 4;
switch cube
case 1
    i = nelx/2;
    j = 1;
    k = 1;
case 2
    i = nelx/2;
    j = 2;
    k = 15;        
case 3
    i = nelx;
    j = 9;14;
    k = nelz;    
case 4    
    i = nelx-2;
    j = 1;
    k = nelz;
end

[~, ~, edofMat_local] = prepare_assemble(microx-1, microx-1, microx-1);

step = microx;
for cube = 1
    switch cube
    case 1
        i = 23;20;
        j = 1;
        k = 1;
    case 2
        i = nelx/2;
        j = 3;
        k = 15;        
    case 3
        i = nelx/2;
        j = 2;
        k = nelz;
    case 4    
        i = nelx-2;
        j = 1;
        k = nelz;
    end
    
    xi = x((j-1)*step+1 : j*step-1, (i-1)*step+1 : i*step-1, (k-1)*step+1 : k*step-1);
    coorx = defx((j-1)*step+1 : j*step, (i-1)*step+1 : i*step, (k-1)*step+1 : k*step);
    coory = defy((j-1)*step+1 : j*step, (i-1)*step+1 : i*step, (k-1)*step+1 : k*step);
    coorz = defz((j-1)*step+1 : j*step, (i-1)*step+1 : i*step, (k-1)*step+1 : k*step);
    def_local= [coorx(:), coory(:), coorz(:)];

    for mat = 1:4
        idx = find(xi==mat);
        if ~isempty(idx)
            [~, def_nodes_new, eles] = separate_rho(xi, edofMat_local, def_local, def_local, mat);
            data_path = ['data/point-clouds/def/local_cube/', num2str(i),'_', num2str(j),'_', num2str(k),'_', num2str(mat),'.vtk'];
            voxel2vtk(data_path, def_nodes_new, eles);        
        end

        % idx = find(xi==mat);
        % if ~isempty(idx)
        %     cx = coorx(idx);
        %     cy = coory(idx);
        %     cz = coorz(idx);
        %     coor = [cx(:), cy(:), cz(:)];
        %     % figure; scatter3(coor(:,1), coor(:,2), coor(:,3));

        %     data_path = ['data/point-clouds/def/def_', 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