三维拓扑优化169行代码

% AN 169 LINE 3D TOPOLOGY OPITMIZATION CODE BY LIU AND TOVAR （JUL 2013）
function [x U c dc]=top3d169（nelxnelynelzvolfracpenalrminmaxloop）
% USER-DEFINED LOOP PARAMETERS
% maxloop = 20;    % Maximum number of iterations
tolx = 0.01;      % Terminarion criterion
displayflag = 0;  % Display structure flag
% USER-DEFINED MATERIAL PROPERTIES
E0 = 1;           % Young‘s modulus of solid material
Emin = 1e-9;      % Young‘s modulus of void-like material
nu = 0.3;         % Poisson‘s ratio
% USER-DEFINED LOAD DOFs
il = nelx; jl = 0; kl = 0:nelz;                         % Coordinates
loadnid = kl*（nelx+1）*（nely+1）+il*（nely+1）+（nely+1-jl）; % Node IDs
loaddof = 3*loadnid（:） - 1;                             % DOFs
% USER-DEFINED SUPPORT FIXED DOFs
[jfkf] = meshgrid（1:nely+11:nelz+1）;                  % Coordinates
fixednid = （kf-1）*（nely+1）*（nelx+1）+jf;                 % Node IDs
fixeddof = [3*fixednid（:）; 3*fixednid（:）-1; 3*fixednid（:）-2]; % DOFs
% PREPARE FINITE ELEMENT ANALYSIS
nele = nelx*nely*nelz;
ndof = 3*（nelx+1）*（nely+1）*（nelz+1）;
F = sparse（loaddof1-1ndof1）;
U = zeros（ndof1）;
freedofs = setdiff（1:ndoffixeddof）;
KE = lk_H8（nu）;
nodegrd = reshape（1:（nely+1）*（nelx+1）nely+1nelx+1）;
nodeids = reshape（nodegrd（1:end-11:end-1）nely*nelx1）;
nodeidz = 0:（nely+1）*（nelx+1）:（nelz-1）*（nely+1）*（nelx+1）;
nodeids = repmat（nodeidssize（nodeidz））+repmat（nodeidzsize（nodeids））;
edofVec = 3*nodeids（:）+1;
edofMat = repmat（edofVec124）+ ...
    repmat（[0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1 ...
    3*（nely+1）*（nelx+1）+[0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1]]nele1）;
iK = reshape（kron（edofMatones（241））‘24*24*nele1）;
jK = reshape（kron（edofMatones（124））‘24*24*nele1）;
% PREPARE FILTER
iH = ones（nele*（2*（ceil（rmin）-1）+1）^21）;
jH = ones（size（iH））;
sH = zeros（size（iH））;
k = 0;
for k1 = 1:nelz
    for i1 = 1:nelx
        for j1 = 1:nely
            e1 = （k1-1）*nelx*nely + （i1-1）*nely+j1;
            for k2 = max（k1-（ceil（rmin）-1）1）:min（k1+（ceil（rmin）-1）nelz）
                for i2 = max（i1-（ceil（rmin）-1）1）:min（i1+（ceil（rmin）-1）nelx）
                    for j2 = max（j1-（ceil（rmin）-1）1）:min（j1+（ceil（rmin）-1）nely）
                        e2 = （k2-1）*nelx*nely + （i2-1）*nely+j2;
                        k = k+1;
                        iH（k） = e1;
                        jH（k） = e2;
                        sH（k） = max（0rmin-sqrt（（i1-i2）^2+（j1-j2）^2+（k1-k2）^2））;
                    end
                end
            end
        end
    end
end
H = sparse（iHjHsH）;
Hs = sum（H2）;
% INITIALIZE ITERATION
x = repmat（volfrac[nelynelxnelz]）;
xPhys = x; 
loop = 0; 
change = 1;
% START ITERATION
while change > tolx && loop < maxloop
    loop = loop+1;
    % FE-ANALYSIS
    sK = reshape（KE（:）*（Emin+xPhys（:）‘.^penal*（E0-Emin））24*24*nele1）;
    K = sparse（iKjKsK）; K = （K+K‘）/2;
    U（freedofs:） = K（freedofsfreedofs）\F（freedofs:）;
    % objectIVE FUNCTION AND SENSITIVITY ANALYSIS
    ce = reshape（sum（（U（edofMat）*KE）.*U（edofMat）2）[nelynelxnelz]）