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.

54 lines
1.5 KiB

3 years ago
%--------------------------
% @Author: Jingqiao Hu
% @Date: 2021-12-07 09:13:42
% @LastEditTime: 2022-01-30 10:12:50
%--------------------------
function [dc, dv] = sensitivity_cvt(c_dH, c_interEles, c_vol, c_rho, c_psi, c_sK, ...
c_edofMat, U, e0, e1, penal, vol0)
ncvt = size(c_rho, 1);
dc_local = cell(ncvt, 1);
parfor ele = 1 : ncvt
%% compliance sensitivity
psi1 = c_psi{ele}; % [12, 3v, m]
psi2 = permute(psi1, [2,1,3]); % [3v, 12, m]
nvdofs = size(psi2, 1);
sK0 = c_sK{ele}; % [12^2, m]
rho = c_rho{ele};
sK1 = sK0 .* (penal * rho(:)'.^(penal-1) * (e1 - e0));
sK2 = reshape(sK1, 12,12,[]); % [12,12,m]
Ke0 = mtimesx(mtimesx(psi2, sK2), psi1); % [3v,3v,m]
Ke1 = reshape(Ke0, nvdofs, []); % [3v, 3v*m]
edof = c_edofMat{ele};
ue = U(edof); % [3v,1]
dc1 = ue' * Ke1; % [1, 3v*m]
dc2 = reshape(dc1, nvdofs, [])'; % [m,3v]
dc_local{ele} = dc2 * ue; % [m, 1]
end
dcdp = zeros(ncvt, 4);
dvdp = zeros(ncvt, 4);
parfor i = 1 : ncvt
for j = 1:4
dH_ele = c_dH{i,j};
interEles = c_interEles{i,j};
nnb = length(interEles);
for k = 1 : nnb
ele = interEles(k);
dcdp(i, j) = dcdp(i, j) + sum(dc_local{ele} .* dH_ele{k});
dvdp(i, j) = dvdp(i, j) + sum(c_vol{ele} / vol0 .* dH_ele{k});
end
end
end
dc = reshape(dcdp', [], 1);
dv = reshape(dvdp', [], 1);
end