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.
65 lines
1.9 KiB
65 lines
1.9 KiB
%--------------------------
|
|
% @Author: Jingqiao Hu
|
|
% @Date: 2021-12-07 09:13:42
|
|
% @LastEditTime: 2021-12-12 12:36:07
|
|
%--------------------------
|
|
function [dc, dv] = sensitivity_cvt_mmc(dH_cell, area_cell, rho_cell, psi_cell, sK_cell, ...
|
|
v_local, edofMat, U, e0, e1, penal, optVcons, area0, nc_local, nd)
|
|
|
|
ncvt = size(rho_cell, 1);
|
|
|
|
pn = 16; % p-norm for volume fraction
|
|
dpnorm = 1/ncvt * (mean(v_local.^pn))^(1/pn-1) * (v_local.^(pn-1)); % [N, 1]
|
|
|
|
nc_global = sum(nc_local) * nd;
|
|
dc = zeros(nc_global, 1);
|
|
dv = zeros(nc_global, 1);
|
|
|
|
% for local fine mesh
|
|
nc_local(2:ncvt+1) = nc_local;
|
|
nc_local(1) = 0;
|
|
|
|
for ele = 1 : ncvt
|
|
%% compliance sensitivity
|
|
psi1 = psi_cell{ele}; % [6, 2v, m]
|
|
psi2 = permute(psi1, [2,1,3]); % [2v, 6, m]
|
|
nvdofs = size(psi2, 1);
|
|
|
|
sK0 = sK_cell{ele}; % [36, m]
|
|
rho = rho_cell{ele};
|
|
|
|
sK1 = sK0 .* (penal * rho(:)'.^(penal-1) * (e1 - e0));
|
|
sK2 = reshape(sK1, 6,6,[]); % [6,6,m]
|
|
|
|
Ke0 = mtimesx(mtimesx(psi2, sK2), psi1); % [2v,2v,m]
|
|
Ke1 = reshape(Ke0, nvdofs, []); % [2v, 2v*m]
|
|
|
|
edof = edofMat{ele};
|
|
ue = U(edof); % [2v,1]
|
|
|
|
dc1 = ue' * Ke1; % [1, 2v*m]
|
|
dc2 = reshape(dc1, nvdofs, [])'; % [m,2v]
|
|
dc3 = dc2 * ue; % [m, 1]
|
|
|
|
m = length(dc3);
|
|
dH0 = dH_cell{ele}; % [nc, 3, m]
|
|
dH1 = permute(dH0, [2,1,3]); % [3, nc, m]
|
|
dH2 = reshape(dH1, [], m); % [3*nc, m]
|
|
|
|
dc4 = dH2 * dc3; % dc_dmmc, [3*nc]
|
|
|
|
idx_begin = sum(nc_local(1:ele))*nd + 1;
|
|
idx_end = sum(nc_local(1:ele+1))*nd;
|
|
dc(idx_begin : idx_end) = dc4;
|
|
|
|
%% volume
|
|
area = area_cell{ele};
|
|
switch optVcons
|
|
case 'pnorm'
|
|
dv0 = dpnorm(ele) * area / sum(area);
|
|
case 'overall'
|
|
dv0 = area / area0;
|
|
end
|
|
dv(idx_begin : idx_end) = dH2 * dv0(:);
|
|
end
|
|
end
|