%-------------------------- % @Author: Jingqiao Hu % @Date: 2021-08-25 17:26:51 % @LastEditTime: 2021-09-09 14:25:42 %-------------------------- function [f0val, df0dx, fval, dfdx] = obj_cons_sens_micro(nelx, nely, microx, H_global, Hs_global, fext, U, nele_global, v0_micro, ... dk_micro, edofMat_ma, rho_cell, cell2global_id, e0, e1, penal, ke1, iK_mi, jK_mi, edofMat_mi, bdofs_num, dD, dofid, newdofs, ... bdofs_plus, ib_dofs, ii_dofs, regulated_matrix, bezier_B, dFdx, De) pnorm = 16; nele_macro = nelx * nely; nele_micro = microx^2; %% compliance & volume_macro comp = fext' * U; % compliance_sensitivities(nelx, nely, microx, dk_micro, U, edofMat_ma, 1, 1); % pnorm-volume constraint of micro volume_micro = reshape(mean(rho_cell, 1), [], 1); % [m,M] -> [1,M] cons_v_micro = (mean(volume_micro.^pnorm))^(1/pnorm) / v0_micro - 1; dpnorm = 1/nele_macro * (mean(volume_micro.^pnorm))^(1/pnorm-1) * (volume_micro.^(pnorm-1)) / v0_micro; dpn_expand = kron(ones(nele_micro, 1), dpnorm); % expand to [nele_global, 1] dv_micro = ones(nele_global, 1) / nele_global .* dpn_expand; % cons_v_micro = mean(rho_cell(:)) / v0_micro - 1; % dv_micro = ones(nele_global, 1) / nele_global / v0_micro; [dc_micro] = obj_sens_micro(nelx, nely, microx, dk_micro, U, edofMat_ma, cell2global_id); dc_micro1 = compliance_sensitivities(nelx, nely, microx, regulated_matrix, U, edofMat_ma, edofMat_mi, dofid, newdofs, ... e0, e1, penal, ke1, iK_mi, jK_mi, bdofs_num, dD, rho_cell, bdofs_plus, ib_dofs, ii_dofs, cell2global_id, bezier_B, dFdx, De); max(dc_micro(:) - dc_micro1(:)) dc_micro(:) = H_global * (dc_micro(:) ./ Hs_global); dv_micro(:) = H_global * (dv_micro(:) ./ Hs_global ); % CHECK: maybe need to be scaled to 1000 times %% sensitivities f0val = comp; df0dx = - dc_micro(:); % note the negetive sign fval = cons_v_micro; dfdx = dv_micro(:)'; end % newdofs = [boundarys, inner]; function dc_micro = compliance_sensitivities(nelx, nely, microx, regulated_matrix, U, edofMat_ma, edofMat_mi, dofid, newdofs, ... e0, e1, penal, ke1, iK_mi, jK_mi, bdofs_num, dD, rho, bdofs_plus, ib_dofs, ii_dofs, cell2global_id, bezier_B, dFdx, De) nele_macro = nelx * nely; nele_micro = microx^2; UM = U(edofMat_ma); % [M,v] order1 = dofid(edofMat_mi'); % order as [border,idofs] order2 = node2ele_id(dofid(edofMat_mi)); % O2 = zeros(2*(microx+1)^2 - bdofs_num, 1); O = zeros(bdofs_num, 1); idofs_plus = 1 - bdofs_plus; % [8,m] dc = zeros(nele_micro, nele_macro); dc1 = zeros(nele_micro, nele_macro); for ele = 1: nele_macro rho_e = rho(:, ele); rho_penal = e0 + rho_e.^penal * (e1-e0); % m,1 micro-rho need to be penal drho = penal * rho_e.^(penal-1) * (e1 - e0); ue = UM(ele, :); % [1,v] Re = regulated_matrix{ele}; % [2n,v] former1 = Re * ue'; % [2n,1] former3 = former1(order1); % expand to [8,m] former4 = former3' * dD .* rho_penal; % [m,8] former5 = reshape(former4', [], 1); former6 = former5(order2); % assemble to [2n,1] former_sub = former6(bdofs_num + 1 : end); % choose dofs corresponding to k_22 of [2i,1] K = assemble_micro_k(dofid, rho_e, microx, e0, e1, penal, ke1, iK_mi, jK_mi); k22 = K(bdofs_num + 1 : end, bdofs_num + 1 : end); % [2i,2i] former_sub1 = k22 \ former_sub; % [2i,1] former7 = [O; former_sub1]'; % [1,2n] former = former7(order1); % [8,m] latter_b(:,1,:) = former3 .* bdofs_plus; % [8,1,m] latter_i(:,1,:) = - former3 .* idofs_plus; % [8,1,m] dke = ke1(:) * drho(:)'; % [8*8, m] dke1 = reshape(dke, 8, 8, nele_micro); % [8, 8, m] dk21 = dke1 .* ib_dofs; % [8, 8, m] dk22 = dke1 .* ii_dofs; % [8, 8, m] dM = mtimesx(dk22, latter_i) - mtimesx(dk21, latter_b); % [8,1,m] dM1 = squeeze(dM); % [8,m] % dM1 = reshape(dM1, [], 1); % dM2 = dM1(order2); % assemble to [2n,1] % dM2_sub = dM2(bdofs_num + 1 : end); % [2i,1] % latter_sub = k22 \ dM2_sub; % [2i,1] % % latter = [O; latter_sub]; % latter = latter(order1); % [8,m] dc(:,ele) = 2 * sum(former .* dM1, 1); % [1,m] % dc(:,ele) = 2 * sum(former4' .* latter, 1); % [1,m] %% ndofs = 2 * (microx + 1)^2; dk = zeros(ndofs, ndofs); O = zeros(bdofs_num, 1); for ele_mi = 12:nele_micro edof = dofid(edofMat_mi(ele_mi, :)); dke = ke1 * drho(ele_mi); dk(edof, edof) = dke; % assemble by the order of [dofid] dk21 = dk(bdofs_num+1 : end, 1 : bdofs_num) * former1(1:bdofs_num); dk22 = dk(bdofs_num+1 : end, bdofs_num+1 : end) * (-former1(bdofs_num+1:end)); dk2_pre = k22 * (dk22 - dk21); % NOTE: not including \bO dphi = [O; dk2_pre]; % a = dphi(edof); % % test = ke1 \ dM1(:,ele_mi); % a = former1(edof)' * dD * rho_penal(ele_mi); % max(dM1(:,ele_mi) - a) dc1(ele_mi,ele) = 2 * former1(edof)' * dD * rho_penal(ele_mi) * dphi(edof); % + former1(edof)' * dD * drho(ele_mi) * former1(edof); % dc2(ele_mi,ele) = 2 * former6(edof)' * dphi(edof); dk(:) = 0; end end dc_micro = dc(cell2global_id); dc_micro1 = dc1(cell2global_id); end function order_id = node2ele_id(edofMat_mi) order1 = reshape(edofMat_mi', [], 1); order1(:,2) = 1:size(order1,1); order2 = sortrows(order1, 1); [~,ia,ic] = unique(order2(:,1)); % C = A(ia) idx = order2(:,2); order_id = idx(ia); end % compliance_micro sensitivities function [dc_micro] = obj_sens_micro(nelx, nely, microx, dk_micro, U, edofMat_ma, cell2global_id) nele_micro = microx^2; nele_macro = nelx * nely; vdofs_num = size(edofMat_ma, 2); % tic UM = U(edofMat_ma); % [M,v] UM1 = kron(ones(nele_micro,1), UM)'; % [v,m*M] UM2 = reshape(UM1, vdofs_num, 1, nele_macro, nele_micro); % [v,1,M,m] UM3 = permute(UM2, [1,2,4,3]); % [v,1,m,M] UM4 = permute(UM3, [2,1,3,4]); % [1,v,m,M] dc5 = squeeze(mtimesx(mtimesx(UM4, dk_micro), UM3)); % [m,M] dc_micro = dc5(cell2global_id); % toc % tic % dc = zeros(nele_micro, nele_macro); % parfor ele = 1: nelx * nely % ue = UM(ele, :); % % dke = dk_micro(:,:,:,ele); % [v,v,m] % dk1 = reshape(dke, vdofs_num, []); % [v,v*m] % % dc_former = ue * dk1; % % [1,v*m] % dc1 = reshape(dc_former, vdofs_num, nele_micro); % % dc(:,ele) = dc1' * ue'; % end % dc_micro1 = dc(cell2global_id); % toc % % max(dc_micro(:) - dc_micro1(:)) end