%?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?????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]）