99行拓扑优化代码学习，转载于另外两篇博客

‘’‘Matlab
%%%% A 99 LINE TOPOLOGY OPTIMIZATION CODE BY OLE SIGMUND, JANUARY 2000 %%%
function top99_notation(nelx,nely,volfrac,penal,rmin)
% INITIALIZE
x(1:nely,1:nelx) = volfrac;
dc(1:nely,1:nelx) = 0.0;
loop = 0; 
change = 1.;
% start ineration 
while change>0.01 
    loop=loop+1;  % loop是迭代计数
    xold=x;       % 将设计变量保存在 变量xold中
    % FE analysis 
    [U]=FE(nelx,nely,x,penal);  % 有限元模型分析，将计算得到的各个节点位移值保存在数组U中 
    % objective function and sensitivity analysis 
    [KE]=lk;   % 计算单元刚度阵 
    c=0.;      % 保存目标函数(柔度)的变量，
    %  遍历所有单元
    % 总的单元数量：nely * nelx;
    % 网格划分数量：（nely+1）*(nelx+1)
    for ely=1:nely 
        for elx=1:nelx 
            % 根据单元在设计域中位置 计算单元节点编号
            n1=(nely+1)*(elx-1)+ely;  % 左上角节点编号  
            n2=(nely+1)*elx    +ely;  % 右上角节点编号  
            % 根据以上两个节点的编号 推演出 单元所有节点的自由度编号
            % 这就限制了此程序只能使用四边形单元剖分四边设计域所得到的有限元模型
            %
            % 从数组 U 中根据节点自由度编号 提取该单元的 位移向量 
            % 左上，右上，右下，左下 自由度
            % 一个点有两个，所以要*2。第一个从1开始，所以*2之后要-1。
            %所示单元的自由度，左上，右上，右下，左下，注意顺序
            Ue=U([2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2],1);
            %  计算单元柔度，叠加到目标函数变量中 
            %  使用SIMP(Solid Isotropic Material with Penalization)材料插值模型
            c=c+x(ely,elx)^penal*Ue'*KE*Ue; 
            %  计算灵敏度  
            dc(ely,elx)=-penal*x(ely,elx)^(penal-1)*Ue'*KE*Ue; 
        end 
    end 
    % filtering of sensitivities 
    [dc]=check(nelx,nely,rmin,x,dc); 
    % design update by the optimality criteria method 
    [x]=OC(nelx,nely,x,volfrac,dc); 
    % print result 
    change = max(max(x-xold));
    disp([' It.: ' sprintf('%4i',loop) ' Obj.: ' sprintf('%10.4f',c) ...
           ' Vol.: ' sprintf('%6.3f',sum(sum(x))/(nelx*nely)) ...
            ' ch.: ' sprintf('%6.3f',change )]);
    % plot densities 
    colormap(gray);
    imagesc(-x);
    axis equal;
    axis tight; 
    axis off;
    pause(1e-6); 
end 
%%%%%%%%%% OPTIMALITY CRITERIA UPDATE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xnew]=OC(nelx,nely,x,volfrac,dc)  
l1 = 0; l2 = 100000; move = 0.2;
while (l2-l1 > 1e-4)
  lmid = 0.5*(l2+l1);
  
  xnew = max(0.001,max(x-move,min(1.,min(x+move,x.*sqrt(-dc./lmid)))));
  if sum(sum(xnew)) - volfrac*nelx*nely > 0
    l1 = lmid;
  else
    l2 = lmid;
  end
end
%%%%%%%%%% MESH-INDEPENDENCY FILTER %%%%%%%%%%%%%%%%%%%
function [dcn]=check(nelx,nely,rmin,x,dc)
% dcn 清零，dcn 用来保存 更新的 目标函数灵敏度
dcn=zeros(nely,nelx);
% 遍历所有单元, 加权平均的一个过程
for i = 1:nelx
  for j = 1:nely
    sum=0.0;
    % 遍历于这个单元相邻的单元
    for k = max(i-floor(rmin),1):min(i+floor(rmin),nelx)
      for l = max(j-floor(rmin),1):min(j+floor(rmin),nely)
        %  sqrt((i-k)^2+(j-l)^2) 是计算此单元与相邻单元的距离
        
        fac = rmin-sqrt((i-k)^2+(j-l)^2);
        sum = sum + max(0,fac);

        dcn(j,i) = dcn(j,i) + max(0,fac)*x(l,k)*dc(l,k);
      end
    end

    dcn(j,i) = dcn(j,i)/(x(j,i)*sum);
  end
end

%%%%%%%%%% FE-ANALYSIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FE analysis 
function [U]=FE(nelx,nely,x,penal) %自定义函数，最后返回[U]
[KE]=lk; %计算单元刚度矩阵  
% sparse 把一个全矩阵转化为一个稀疏矩阵，只存储每一个非零元素的三个值：元素值，元素的行号和列号
% 总体刚度矩阵的稀疏矩阵
% *2是因为x，y都有一个数
% F = KU
K=sparse(2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1)); 
F=sparse(2*(nely+1)*(nelx+1),1); % 外力
U=sparse(2*(nely+1)*(nelx+1),1); % 要求的结果u
%组装整体刚度矩阵  
for elx=1:nelx 
    for ely=1:nely 
        %一列列的排序，y*x
        n1=(nely+1)*(elx-1)+ely; % 左上
        n2=(nely+1)*elx    +ely; % 右上
        % 左上，右上，右下，左下 自由度
        % 一个点有两个，所以要*2。第一个从1开始，所以*2之后要-1。
        edof=[2*n1-1;2*n1;2*n2-1;2*n2;2*n2+1;2*n2+2;2*n1+1;2*n1+2]; 
        %将单元刚度矩阵组装成总的刚度矩阵
        K(edof,edof)=K(edof,edof)+x(ely,elx)^penal*KE; 
    end 
end 
% define loads and supports 
 ip=(nelx+1)*(nely+1); 
 F(2*ip,1)=-1;  % 施加载荷 % 应用了一个在右下角的垂直单元力。             

% 施加位移约束  
fixeddofs = (1:2*(nely+1)); 
alldofs   = (1:2*(nely+1)*(nelx+1)); 
% 所有不受约束的节点自由度 
freedofs  =setdiff(alldofs,fixeddofs); 
% solving 求解线性方程组, 得到节点自由度的位移值
% 注意在此仅对 freedofs做计算，原理见！！！那里
U(freedofs,:)=K(freedofs,freedofs)\F(freedofs,:); 
U(fixeddofs,:)=0;  % 受约束的节点自由度的位移值 设为  0 

%%%%%%%%%% ELEMENT STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [KE]=lk
E = 1.; % 杨氏模量，直观可理解为这个材料的硬度，越大则越硬
nu = 0.3; % 泊松比，当这个材料受到压缩时，nu>0则会膨胀（常见），nu<0则会收缩
% 预计算了一个矩形的单元刚度矩阵，具体求解可以看曾攀的《有限元分析》
k=[ 1/2-nu/6   1/8+nu/8 -1/4-nu/12 -1/8+3*nu/8 ... 
   -1/4+nu/12 -1/8-nu/8  nu/6       1/8-3*nu/8];
KE = E/(1-nu^2)*[ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8)
                  k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3)
                  k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2)
                  k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5)
                  k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4)
                  k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7)
                  k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6)
                  k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1)];
              

原始的博文地址：

Sigmund的99行Matlab拓扑优化程序简析_cocoonyang的专栏-CSDN博客_99行拓扑优化

[更新] 99行拓扑优化 代码解析_我要上火星（停止更新，转为私人笔记）-CSDN博客_拓扑优化代码