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.
42 lines
1.2 KiB
42 lines
1.2 KiB
3 years ago
|
%--------------------------
|
||
|
|
% @Author: Jingqiao Hu
|
||
|
|
% @Date: 2021-11-28 13:49:21
|
||
|
|
% @LastEditTime: 2022-02-01 12:02:45
|
||
|
|
%--------------------------
|
||
|
|
function K = stiffness_cvt(rho_cell, sK_cell, psi_cell, c_edofMat, nvar, e0, e1, penal, c_vol)
|
||
|
|
|
||
|
|
ncvt = size(psi_cell, 1);
|
||
|
|
K = sparse(nvar, nvar);
|
||
|
|
|
||
|
|
for ele = 1 : ncvt
|
||
|
|
psi1 = psi_cell{ele}; % [12, 3v, m]
|
||
|
|
psi2 = permute(psi1, [2,1,3]); % [3v, 12, m]
|
||
|
|
|
||
|
|
sK0 = sK_cell{ele}; % [12^2, m]
|
||
|
|
rho = rho_cell{ele};
|
||
|
|
|
||
|
|
sK1 = sK0 .* (e0 + rho(:)'.^penal * (e1 - e0));
|
||
|
|
sK2 = reshape(sK1, 12,12,[]); % [12,12,m]
|
||
|
|
|
||
|
|
Ke0 = mtimesx(mtimesx(psi2, sK2), psi1); % [3v,3v,m]
|
||
|
|
|
||
|
|
% vol0 = c_vol{ele}; % [m,1]
|
||
|
|
% vol1 = zeros(1,1,length(vol0));
|
||
|
|
% vol1(1,1,:) = vol0;
|
||
|
|
% Ke0 = Ke0 .* vol1;
|
||
|
|
|
||
|
|
Ke1 = squeeze(sum(Ke0, 3)); % [3v,3v], each ele is different
|
||
|
|
|
||
|
|
edof = c_edofMat{ele};
|
||
|
|
K(edof, edof) = K(edof, edof) + Ke1;
|
||
|
|
end
|
||
|
|
|
||
|
|
% add small value to diagonal elements
|
||
|
|
% in case of disconnected projected rho caused by thin components
|
||
|
|
S = spdiags(repmat(1e-5, nvar, 1),0,nvar,nvar);
|
||
|
|
K = K + S;
|
||
|
|
|
||
|
|
if sum(isinf(K(:))) > 0 || sum(isnan(K(:))) > 0
|
||
|
|
error("here");
|
||
|
|
end
|
||
|
|
end
|