%-------------------------- % @Author: Jingqiao Hu % @Date: 2021-11-02 15:18:36 % @LastEditTime: 2021-11-04 15:56:57 % Ke is [6,6,faces_num] %-------------------------- function optimize_t_trimesh(Ke, nodes, faces, nodes_cell, faces_cell, MMCs, MMCs_max, MMCs_min, radius, p, ... edofMat, iK, jK, freedofs, F, nidx_cell, fidx_cell, MMCs_num, volfrac) nele_global = size(faces, 1); cvt_num = size(MMCs, 1); pn = 16; % for pnorm U = zeros(size(nodes,1)*2, 1); [xval, xold1, xold2, xmin, xmax, low, upp, m_mma, c_mma, d_mma, a0_mma, a_mma] = MMA_paras1(MMCs, MMCs_min, MMCs_max, MMCs_num); loop = 1; change = 1; maxloop = 100; while loop < maxloop && change > 1e-4 %% MMC-part [H, dH_global, denk_global, den_global, v_local] = project_MMC(nodes_cell, faces_cell, MMCs, radius, p, nidx_cell, fidx_cell); figure(1); clf; trisurf(faces, nodes(:,1), nodes(:,2), H); %% simulation [dcdH, comp] = simulation_MMC_trimesh(edofMat, iK, jK, freedofs, U, F, Ke, denk_global, den_global); cons_v_micro = (mean(v_local.^pn))^(1/pn); %% Sensitivities dpnorm = 1/cvt_num * (mean(v_local.^pn))^(1/pn-1) * (v_local.^(pn-1)); % [cvy_num,1] dpn_expand = zeros(nele_global, 1); for i = 1 : cvt_num dpn_expand(fidx_cell(i,1) : fidx_cell(i,2)) = dpnorm(i); end dvdH = ones(nele_global, 1) / nele_global .* dpn_expand; dc = zeros(MMCs_num, 3); dv = zeros(MMCs_num, 3); idx = 0; for i = 1 : cvt_num % NOTE: each MMC only effects its own cell dH = dH_global{i}; edges_num = size(MMCs{i}, 2); dcdt = zeros(edges_num, 3); dvdt = zeros(edges_num, 3); for e = 1 : edges_num dHdt = squeeze(dH(e, :, :))'; % [m,3] dcdt(e, :) = dcdH(:)' * dHdt; % [1,m]*[m,3] -> [1,3] dvdt(e, :) = dvdH(:)' * dHdt; end dc(idx+1 : idx+edges_num, :) = dcdt; dv(idx+1 : idx+edges_num, :) = dvdt; end dc = dc'; % change to [3, MMCs_num] dv = dv'; %% optimization-part f0val= comp; df0dx= -dc(:) / max(abs(dc(:))); fval = cons_v_micro - volfrac; dfdx = dv(:) / max(abs(dv(:))); [xmma, ~, ~, ~, ~, ~, ~, ~, ~,low, upp]= mmasub_mmc(m_mma, nn, Loop, xval(:), xmin, xmax, xold1,... xold2, f0val, df0dx, fval, dfdx', low, upp, a0_mma, a_mma, c_mma, d_mma); xold2 = xold1; xold1 = xval; change = max(abs(xval-xmma)); % update MMCs MMCs = update_MMC(xval, MMCs); end end %% utils function [dcdH, comp] = simulation_MMC_trimesh(edofMat, iK, jK, freedofs, U, F, Ke, denk, den) sKe = reshape(Ke, 36, []); % [36,m] sK = sKe .* denk(:)'; K = sparse(iK(:),jK(:),sK(:)); K=(K+K')/2; U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:); % Energy of element u0 = U(edofMat)'; % [6,m] u1(1,:,:) = u0; % [1,6,m] u2(:,1,:) = u0; % [6,1,m] energy = squeeze(mtimesx(mtimesx(u1, Ke), u2)); % [m,1] dcdH = 2 * energy(:) .* den(:); % [m,1] comp = F'*U; end function [H_global, dH_global, denk_global, den_global, v_local] = project_MMC(nodes_cell, faces_cell, MMCs, radius, p, nidx_cell, fidx_cell) cvt_num = size(MMCs, 1); denk_global = zeros(fidx_cell(end), 1); den_global = zeros(fidx_cell(end), 1); v_local = zeros(cvt_num, 1); % for volume-pnorm H_global = zeros(nidx_cell(end), 1); dH_global = cell(cvt_num, 1); step = radius / 2; for i = 1:cvt_num mmc = MMCs{i}; % [7, edges_num] nodes = nodes_cell{i}; faces = faces_cell{i}; edges_num = size(mmc,2); % edges_num each_phi = cell(edges_num, 1); phi0 = 0; for j = 1 : edges_num tmpPhi = tPhi(mmc(:,j), nodes(:,1), nodes(:,2)); each_phi{j} = tmpPhi; phi0 = phi0 + exp(p * tmpPhi); end phi_max = log(phi0) / p; H = Heaviside_simply(phi_max, radius); H_global(nidx_cell(i,1) : nidx_cell(i,2)) = H; edofMat0 = [2*faces-1, 2*faces]; % [n,6] edofMat = edofMat0(:, [1,4,2,5,3,6]); EleNodesID = edofMat(:, 2:2:6) ./ 2; denk= sum(H(EleNodesID).^2, 2) / 3; % 3 for tiangle den = sum(H(EleNodesID), 2) / 3; % [m,1] denk_global(fidx_cell(i,1) : fidx_cell(i,2)) = denk; den_global(fidx_cell(i,1) : fidx_cell(i,2)) = den; v_local(i) = mean(den); % figure(3); clf; trimesh(faces, nodes(:,1), nodes(:,2), H); faces_num = size(faces, 1); dH_global{i} = sensitivity_MMC(edges_num, mmc, step, nodes, p, each_phi, radius, EleNodesID, faces_num); end end function dH = sensitivity_MMC(edges_num, mmc, step, nodes, p, each_phi, radius, EleNodesID, faces_num) dH = zeros(edges_num, 3, faces_num); for j = 1 : edges_num mmc_j = mmc(:,j); for ti = 1 : 3 mmc_plus = mmc_j; mmc_minus = mmc_j; mmc_plus(ti) = mmc_j(ti) + step; phi_plus0 = tPhi(mmc_plus, nodes(:,1), nodes(:,2)); phi_plus = exp(p * phi_plus0); for j1 = 1 : j-1 phi_plus = phi_plus + exp(p * each_phi{j1}); end for j2 = j+1 : edges_num phi_plus = phi_plus + exp(p * each_phi{j2}); end phi_plus_blend = log(phi_plus) / p; mmc_minus(ti) = mmc_j(ti) - step; phi_minus0 = tPhi(mmc_minus, nodes(:,1), nodes(:,2)); phi_minus = exp(p * phi_minus0); for j1 = 1 : j-1 phi_minus = phi_minus + exp(p * each_phi{j1}); end for j2 = j+1 : edges_num phi_minus = phi_minus + exp(p * each_phi{j2}); end phi_minus_blend = log(phi_minus) / p; H_plus = Heaviside_simply(phi_plus_blend, radius); H_minus = Heaviside_simply(phi_minus_blend, radius); dHi = (H_plus - H_minus) / (2 * step); % [n,1] dH(j, ti, :) = sum(dHi(EleNodesID), 2) / 3; % [n,1] -> [m,1] % CHECK: need dH/dH_e ? from node to ele % figure(1); clf; trimesh(faces, nodes(:,1), nodes(:,2), H_plus); % figure(2); clf; trimesh(faces, nodes(:,1), nodes(:,2), H_minus); end end end function MMCs = update_MMC(xval, MMCs) cvt_num = size(MMCs, 1); idx = 0; for i = 1:cvt_num mmc = MMCs{i}; edges_num = size(mmc, 2); MMCs{i}(1:3, :) = xval(:, idx+1 : idx+edges_num); end end function [xval, xold1, xold2, xmin, xmax, low, upp, m_mma, c_mma, d_mma, a0_mma, a_mma] = MMA_paras1(MMCs, MMCs_min, MMCs_max, MMCs_num) xval = zeros(3, MMCs_num); xmin = zeros(3, MMCs_num); xmax = zeros(3, MMCs_num); cvt_num = size(MMCs, 1); idx = 0; for i = 1 : cvt_num mmc = MMCs{i}; edges_num = size(mmc, 2); xval(:, idx+1 : idx+edges_num) = mmc(1:3, :); xmin(:, idx+1 : idx+edges_num) = MMCs_min{i}(1:3, :); xmax(:, idx+1 : idx+edges_num) = MMCs_max{i}(1:3, :); end % variables are t1,t2,t3 xval = xval(:); xold1= xval; xold2 = xval; xmin = xmin(:); xmax = xmax(:); low = xmin; upp = xmax; m_mma = 1; % number of constrainm c_mma = 1000*ones(m_mma,1); d_mma = zeros(m_mma,1); a0_mma = 1; a_mma = zeros(m_mma,1); end