前言

Federico Ferrari 和Ole Sigmund的高效3D拓扑优化程序，非常值的学习。

图片摘自于《A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D》

一、top3D125

top3D125可以实现高效3D结构拓扑优化

二、使用步骤

top3D125(24,12,12,0.12,3,sqrt(3),1,‘N’,0.5,1,0.2,100);

三、主程序

function top3D125(nelx,nely,nelz,volfrac,penal,rmin,ft,ftBC,eta,beta,move,maxit)
%top3D125(24,12,12,0.12,3,sqrt(3),1,'N',0.5,1,0.2,100);
% ---------------------------- PRE. 1) MATERIAL AND CONTINUATION PARAMETERS
E0 = 1;                                                                    % Young modulus of solid
Emin = 1e-9;                                                               % Young modulus of "void"
nu = 0.3;                                                                  % Poisson ratio
penalCnt = { 1, 1, 25, 0.25 };                                             % continuation scheme on penal
betaCnt  = { 1, 1, 25,    2 };                                             % continuation scheme on beta
if ftBC == 'N', bcF = 'symmetric'; else, bcF = 0; end                      % filter BC selector
% ----------------------------------------- PRE. 2) DISCRETIZATION FEATURES
nEl = nelx * nely * nelz;                                                  % number of elements          #3D#
nodeNrs = int32( reshape( 1 : ( 1 + nelx ) * ( 1 + nely ) * ( 1 + nelz ), ...1 + nely, 1 + nelz, 1 + nelx ) );                                      % nodes numbering             #3D#
cVec = reshape( 3 * nodeNrs( 1 : nely, 1 : nelz, 1 : nelx ) + 1, nEl, 1 ); %                             #3D#
cMat = cVec+int32( [0,1,2,3*(nely+1)*(nelz+1)+[0,1,2,-3,-2,-1],-3,-2,-1,3*(nely+...1)+[0,1,2],3*(nely+1)*(nelz+2)+[0,1,2,-3,-2,-1],3*(nely+1)+[-3,-2,-1]]);% connectivity matrix         #3D#
nDof = ( 1 + nely ) * ( 1 + nelz ) * ( 1 + nelx ) * 3;                     % total number of DOFs        #3D#
[ sI, sII ] = deal( [ ] );
for j = 1 : 24sI = cat( 2, sI, j : 24 );sII = cat( 2, sII, repmat( j, 1, 24 - j + 1 ) );
end
[ iK , jK ] = deal( cMat( :,  sI )', cMat( :, sII )' );
Iar = sort( [ iK( : ), jK( : ) ], 2, 'descend' ); clear iK jK              % reduced assembly indexing
Ke = 1/(1+nu)/(2*nu-1)/144 *( [ -32;-6;-6;8;6;6;10;6;3;-4;-6;-3;-4;-3;-6;10;...3;6;8;3;3;4;-3;-3; -32;-6;-6;-4;-3;6;10;3;6;8;6;-3;-4;-6;-3;4;-3;3;8;3;...3;10;6;-32;-6;-3;-4;-3;-3;4;-3;-6;-4;6;6;8;6;3;10;3;3;8;3;6;10;-32;6;6;...-4;6;3;10;-6;-3;10;-3;-6;-4;3;6;4;3;3;8;-3;-3;-32;-6;-6;8;6;-6;10;3;3;4;...-3;3;-4;-6;-3;10;6;-3;8;3;-32;3;-6;-4;3;-3;4;-6;3;10;-6;6;8;-3;6;10;-3;...3;8;-32;-6;6;8;6;-6;8;3;-3;4;-3;3;-4;-3;6;10;3;-6;-32;6;-6;-4;3;3;8;-3;...3;10;-6;-3;-4;6;-3;4;3;-32;6;3;-4;-3;-3;8;-3;-6;10;-6;-6;8;-6;-3;10;-32;...6;-6;4;3;-3;8;-3;3;10;-3;6;-4;3;-6;-32;6;-3;10;-6;-3;8;-3;3;4;3;3;-4;6;...-32;3;-6;10;3;-3;8;6;-3;10;6;-6;8;-32;-6;6;8;6;-6;10;6;-3;-4;-6;3;-32;6;...-6;-4;3;6;10;-3;6;8;-6;-32;6;3;-4;3;3;4;3;6;-4;-32;6;-6;-4;6;-3;10;-6;3;...-32;6;-6;8;-6;-6;10;-3;-32;-3;6;-4;-3;3;4;-32;-6;-6;8;6;6;-32;-6;-6;-4;...-3;-32;-6;-3;-4;-32;6;6;-32;-6;-32]+nu*[ 48;0;0;0;-24;-24;-12;0;-12;0;...24;0;0;0;24;-12;-12;0;-12;0;0;-12;12;12;48;0;24;0;0;0;-12;-12;-24;0;-24;...0;0;24;12;-12;12;0;-12;0;-12;-12;0;48;24;0;0;12;12;-12;0;24;0;-24;-24;0;...0;-12;-12;0;0;-12;-12;0;-12;48;0;0;0;-24;0;-12;0;12;-12;12;0;0;0;-24;...-12;-12;-12;-12;0;0;48;0;24;0;-24;0;-12;-12;-12;-12;12;0;0;24;12;-12;0;...0;-12;0;48;0;24;0;-12;12;-12;0;-12;-12;24;-24;0;12;0;-12;0;0;-12;48;0;0;...0;-24;24;-12;0;0;-12;12;-12;0;0;-24;-12;-12;0;48;0;24;0;0;0;-12;0;-12;...-12;0;0;0;-24;12;-12;-12;48;-24;0;0;0;0;-12;12;0;-12;24;24;0;0;12;-12;...48;0;0;-12;-12;12;-12;0;0;-12;12;0;0;0;24;48;0;12;-12;0;0;-12;0;-12;-12;...-12;0;0;-24;48;-12;0;-12;0;0;-12;0;12;-12;-24;24;0;48;0;0;0;-24;24;-12;...0;12;0;24;0;48;0;24;0;0;0;-12;12;-24;0;24;48;-24;0;0;-12;-12;-12;0;-24;...0;48;0;0;0;-24;0;-12;0;-12;48;0;24;0;24;0;-12;12;48;0;-24;0;12;-12;-12;...48;0;0;0;-24;-24;48;0;24;0;0;48;24;0;0;48;0;0;48;0;48 ] );             % elemental stiffness matrix  #3D#
Ke0( tril( ones( 24 ) ) == 1 ) = Ke';
Ke0 = reshape( Ke0, 24, 24 );
Ke0 = Ke0 + Ke0' - diag( diag( Ke0 ) );                                    % recover full matrix
% ----------------------------- PRE. 3) LOADS, SUPPORTS AND PASSIVE DOMAINS
lcDof = 3 * nodeNrs( 1 : nely + 1, 1, nelx + 1 );
fixed = 1 : 3 * ( nely + 1 ) * ( nelz + 1 );
[ pasS, pasV ] = deal( [], [] );                                           % passive solid and void elements
F = fsparse( lcDof, 1, -sin((0:nely)/nely*pi)', [ nDof, 1 ] );             % define load vector
free = setdiff( 1 : nDof, fixed );                                         % set of free DOFs
act = setdiff( ( 1 : nEl )', union( pasS, pasV ) );                        % set of active d.v.
% --------------------------------------- PRE. 4) DEFINE IMPLICIT FUNCTIONS
prj = @(v,eta,beta) (tanh(beta*eta)+tanh(beta*(v(:)-eta)))./...(tanh(beta*eta)+tanh(beta*(1-eta)));                                   % projection
deta = @(v,eta,beta) - beta * csch( beta ) .* sech( beta * ( v( : ) - eta ) ).^2 .* ...sinh( v( : ) * beta ) .* sinh( ( 1 - v( : ) ) * beta );                % projection eta-derivative
dprj = @(v,eta,beta) beta*(1-tanh(beta*(v-eta)).^2)./(tanh(beta*eta)+tanh(beta*(1-eta)));% proj. x-derivative
cnt = @(v,vCnt,l) v+(l>=vCnt{1}).*(v<vCnt{2}).*(mod(l,vCnt{3})==0).*vCnt{4};
% -------------------------------------------------- PRE. 5) PREPARE FILTER
[dy,dz,dx]=meshgrid(-ceil(rmin)+1:ceil(rmin)-1,...-ceil(rmin)+1:ceil(rmin)-1,-ceil(rmin)+1:ceil(rmin)-1 );
h = max( 0, rmin - sqrt( dx.^2 + dy.^2 + dz.^2 ) );                        % conv. kernel                #3D#
Hs = imfilter( ones( nely, nelz, nelx ), h, bcF );                         % matrix of weights (filter)  #3D#
dHs = Hs;
% ------------------------ PRE. 6) ALLOCATE AND INITIALIZE OTHER PARAMETERS
[ x, dsK, dV ] = deal( zeros( nEl, 1 ) );                                  % initialize vectors
dV( act, 1 ) = 1/nEl/volfrac;                                              % derivative of volume
x( act ) = ( volfrac*( nEl - length(pasV) ) - length(pasS) )/length( act );% volume fraction on active set
x( pasS ) = 1;                                                             % set x = 1 on pasS set
[ xPhys, xOld, ch, loop, U ] = deal( x, 1, 1, 0, zeros( nDof, 1 ) );       % old x, x change, it. counter, U
% ================================================= START OPTIMIZATION LOOP
while ch > 1e-6 && loop < maxitloop = loop + 1;                                                         % update iter. counter% ----------- RL. 1) COMPUTE PHYSICAL DENSITY FIELD (AND ETA IF PROJECT.)xTilde = imfilter( reshape( x, nely, nelz, nelx ), h, bcF ) ./ Hs;       % filtered field              #3D#xPhys( act ) = xTilde( act );                                            % reshape to column vectorif ft > 1                              % compute optimal eta* with Newtonf = ( mean( prj( xPhys, eta, beta ) ) - volfrac )  * (ft == 3);      % function (volume)while abs( f ) > 1e-6           % Newton process for finding opt. etaeta = eta - f / mean( deta( xPhys, eta, beta ) );f = mean( prj( xPhys, eta, beta ) ) - volfrac;enddHs = Hs ./ reshape( dprj( xPhys, eta, beta ), nely, nelz, nelx );   % sensitivity modification    #3D#xPhys = prj( xPhys, eta, beta );                                     % projected (physical) fieldendch = norm( xPhys - xOld ) ./ nEl;xOld = xPhys;% -------------------------- RL. 2) SETUP AND SOLVE EQUILIBRIUM EQUATIONSsK = ( Emin + xPhys.^penal * ( E0 - Emin ) );dsK( act ) = -penal * ( E0 - Emin ) * xPhys( act ) .^ ( penal - 1 );sK = reshape( Ke( : ) * sK', length( Ke ) * nEl, 1 );K = fsparse( Iar( :, 1 ), Iar( :, 2 ), sK, [ nDof, nDof ] );L = chol( K( free, free ), 'lower' );U( free ) = L' \ ( L \ F( free ) );                                      % f/b substitution% ------------------------------------------ RL. 3) COMPUTE SENSITIVITIESdc = dsK .* sum( ( U( cMat ) * Ke0 ) .* U( cMat ), 2 );                  % derivative of compliancedc = imfilter( reshape( dc, nely, nelz, nelx ) ./ dHs, h, bcF );         % filter objective sens.      #3D#dV0 = imfilter( reshape( dV, nely, nelz, nelx ) ./ dHs, h, bcF );        % filter compliance sens.     #3D#% ----------------- RL. 4) UPDATE DESIGN VARIABLES AND APPLY CONTINUATIONxT = x( act );[ xU, xL ] = deal( xT + move, xT - move );                               % current upper and lower boundocP = xT .* sqrt( - dc( act ) ./ dV0( act ) );                           % constant part in resizing rulel = [ 0, mean( ocP ) / volfrac ];                                        % initial estimate for LMwhile ( l( 2 ) - l( 1 ) ) / ( l( 2 ) + l( 1 ) ) > 1e-4                   % OC resizing rulelmid = 0.5 * ( l( 1 ) + l( 2 ) );x( act ) = max( max( min( min( ocP / lmid, xU ), 1 ), xL ), 0 );if mean( x ) > volfrac, l( 1 ) = lmid; else, l( 2 ) = lmid; endend[penal,beta] = deal(cnt(penal,penalCnt,loop), cnt(beta,betaCnt,loop));   % apply conitnuation on parameters% -------------------------- RL. 5) PRINT CURRENT RESULTS AND PLOT DESIGNfprintf( 'It.:%5i C:%6.5e V:%7.3f ch.:%0.2e penal:%7.2f beta:%7.1f eta:%7.2f lm:%0.2e \n', ...loop, F'*U, mean(xPhys(:)), ch, penal, beta, eta, lmid );isovals = shiftdim( reshape( xPhys, nely, nelz, nelx ), 2 );isovals = smooth3( isovals, 'box', 1 );%patch(isosurface(isovals, .5),'FaceColor','b','EdgeColor','none');%patch(isocaps(isovals, .5),'FaceColor','r','EdgeColor','none');patch(isosurface(isovals, .5),'FaceColor','b','EdgeColor','K');patch(isocaps(isovals, .5),'FaceColor','r','EdgeColor','K');drawnow; view( [ 145, 25 ] ); axis equal tight off; cla();
end
end

四、程序说明

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This Matlab code was written by F. Ferrari, O. Sigmund                   %
% Dept. of Solid Mechanics-Technical University of Denmark,2800 Lyngby (DK)%
% Please send your comments to: feferr@mek.dtu.dk                          %
%                                                                          %
% The code is intended for educational purposes and theoretical details    %
% are discussed in the paper Ferrari, F. Sigmund, O. - A new generation 99 %
% line Matlab code for compliance Topology Optimization and its extension  %
% to 3D, SMO, 2020                                                         %
%                                                                          %
% The code as well as a postscript version of the paper can be             %
% downloaded from the web-site: http://www.topopt.dtu.dk                   %
%                                                                          %
% Disclaimer:                                                              %
% The authors reserves all rights but do not guaranty that the code is     %
% free from errors. Furthermore, we shall not be liable in any event       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

四、备注

程序运行需要安装 [stenglib-master] 包，不然运行不了；另外感谢两位大佬的贡献，从他们那里学到了不少，如有侵权请通知我删帖。
参考文献：《A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D》
参考网址： http://www.topopt.dtu.dk

Federico Ferrari 和Ole Sigmund的高效3D拓扑优化程序相关推荐

matlab拓扑优化流程图,Sigmund的99行Matlab拓扑优化程序简析
引言 Sigmund在2001年在Structural and Multidisciplinary Optimization 上发表一篇名为 "A 99 line topology opti ...
Sigmund的99行Matlab拓扑优化程序简析
引言 Sigmund在2001年在Structural and Multidisciplinary Optimization 发表一篇名为 "A 99 line topology optim ...
python 3d打印_基于Python的结构拓扑优化与3D打印试验研究
收稿日期:2017 -05 -22 修回日期:2017 -06 -01 第 35 卷第 8 期计算机仿真 2018 年 8 月文章编号:1006 -9348( 2018) 08 -017 ...
Arm和Unity联合推出：适用于移动应用程序的3D美术优化-[5]光照
目录 1.光照 - 概述 1.概述 2.关于渲染管线的说明 2.光照模式和类型 1.光照模式 2.烘焙光照 3.实时光照 4.混合光照 5.实时光源和光源类型 6.总结 3.静态对象的光照 1.概述 ...
Arm和Unity联合推出：适用于移动应用程序的3D美术优化-[2]几何体
目录 1.几何体 - 简介 1.1.概述 1.2.什么是几何体? 1.3.总结 2.三角形和多边形的使用量 2.1.概述 2.2.在不同的硬件/屏幕大小上进行测试 2.3.重要区域的细节 2.4.避免 ...
Arm和Unity联合推出：适用于移动应用程序的3D美术优化-[3]纹理
目录 1.纹理-简介 1.1.概述 1.2.什么是纹理 1.3.总结 2.纹理大小.颜色空间和压缩 2.1.概述 2.2.纹理大小 2.3.纹理颜色空间 2.4.纹理压缩 2.5.总结 3.纹理图集 ...
基于HT for Web的3D拓扑树的实现
在HT for Web中2D和3D应用都支持树状结构数据的展示,展现效果各异,2D上的树状结构在展现层级关系明显,但是如果数据量大的话,看起来就没那么直观,找到指定的节点比较困难,而3D上的树状结构在 ...
安全看得见，阿里云性能监控 ARMS 全真3D拓扑实现一“屏”了然
微服务架构下,各类服务之间存在着错综复杂的依赖关系.一旦业务出现问题,追查问题源头就好比大海捞针,没有头绪.但业务不等人,此时,在最短的时间内定位问题根源是开发和运维人员对微服务监控产品的核心诉求. ...
[转]vb高效编程（优化）
本文适合任何水平的vb编程人员. 一.减少加载窗体数目每一个加载的窗体,无论可视与否,都要占据一定数量的内存(其数量随窗体上控件的类型和数量,以及窗体上位图的大小等的不同而变化).只在需要显示时才加 ...

Federico Ferrari 和Ole Sigmund的高效3D拓扑优化程序

文章目录

前言