%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-07-12 09:23:41
 % @LastEditTime: 2021-08-30 19:19: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;

% recover deformation
uglobal = zeros(globaly, globalx, globalz, 3);
uglobal(:, :, :, 1) = defx - gridx;
uglobal(:, :, :, 2) = defy - gridy;
uglobal(:, :, :, 3) = defz - gridz;


step = 10;
for it = 1:8
    
    uglobal_it = uglobal / 10 * it;
    
    defx = uglobal_it(:,:,:,1) + gridx;
    defy = uglobal_it(:,:,:,2) + gridy0;
    defz = uglobal_it(:,:,:,3) + gridz;        

    for i = 1 : (3 * microx / step)
        xi = x((i-1)*step+1 : i*step, :, :);
        coorx = defx((i-1)*step+1 : i*step, :, :);
        coory = defy((i-1)*step+1 : i*step, :, :);
        coorz = defz((i-1)*step+1 : i*step, :, :);

        for mat = 1:4
            idx = find(xi==mat);
            if ~isempty(idx)
                cx = coorx(idx);
                cy = coory(idx);
                cz = coorz(idx);
                
                [f,v] = isosurface(X,Y,Z,V,isovalue);
                
                coor = [cx(:), cy(:), cz(:)];
                % figure; scatter3(coor(:,1), coor(:,2), coor(:,3));

                data_path = ['data/point-clouds/def/iteration/',num2str(it),'/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
end