You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

188 lines
6.8 KiB

3 years ago
%--------------------------
% @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