The code and data of the article " Topology optimization of structures with steady-state heat conduction using an improved parameterized level set method "
1 The table of input parameters and the code of Example 1
|
nelx |
nely |
volfrac |
ep |
dsp |
kmin |
|
60 |
60 |
0.5 |
0.01 |
3 |
0.01 |
function TOCSPRBF_MMA(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=60;volfrac=0.5;ep=0.01;nele = nelx*nely;
%% LEVEL SET FUNCTION INITIALIZATION水平集函数初始化
r=nely*0.05;
hX = nelx*[repmat([1/6,1/3,1/2,2/3,5/6],1,5)];
hY = nely*[kron([1/6,1/3,1/2,2/3,5/6],ones(1,5))];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION径向基函数初始化
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=10000;d=0;m=1;%MMA参数
%% THERMAL FINITE ELEMENT ANALYSIS PREPARATION
TA1 = [ 8 -2; -2 8];
TA2 = [-4 -2; -2 -4];
KEth = 1/12*[TA1 TA2; TA2 TA1];
TeleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
TeleNode = repmat(TeleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
TedofMat = TeleNode;
TiK = reshape(kron(TedofMat,ones(4,1))',16*nele,1);
TjK = reshape(kron(TedofMat,ones(1,4))',16*nele,1);
Tmaxdof = (nely+1)*(nelx+1);
% DEFINE LOADS AND SUPPORTS (HEATED BOTTOM EDGE)
k0 = 1; kmin = 0.01;
T0 = 0;
Q = sparse(((nely+1)*nelx/2+ceil(nely/2)+1),1,10,nNode,1);
T = ones(Tmaxdof,1)*T0;
TBC1 = 1:nely+1:Tmaxdof;
TBC2 = nely+1:nely+1:Tmaxdof;
TBC3 = 1:nely+1;
TBC4 = (nely+1)*nelx+1:(nely+1)*(nelx+1);
Tcdofs = [TBC1 TBC2 TBC3 TBC4];
Tfdofs = setdiff(1:Tmaxdof,Tcdofs);
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 100; nRelax = 10;
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(TeleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(TeleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(TeleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(TeleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
TsK = reshape(KEth(:)*(kmin+eleVol'*(k0-kmin)),16*nelx*nely,1);
Kth = sparse(TiK,TjK,TsK); Kth = (Kth+Kth')/2;
T(Tfdofs) = Kth(Tfdofs,Tfdofs)\(Q(Tfdofs)-Kth(Tfdofs,Tcdofs)*T(Tcdofs));
TeleComp = 1/2*sum((T(TedofMat)*KEth).*T(TedofMat),2).*(kmin+eleVol*(k0-kmin));
comp(iT) = sum(TeleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
TeleComp = reshape(TeleComp,nely,nelx);
TeleCompLR = [TeleComp(:,1),TeleComp]+[TeleComp,TeleComp(:,end)];
TnodeComp = ([TeleCompLR;TeleCompLR(end,:)]+[TeleCompLR(1,:);TeleCompLR])/4;
%% MMA准备
DeltaPhi=reshape(DeltaPhi,nNode,1);
TnodeComp=reshape(TnodeComp,nNode,1);
dJda=-((TnodeComp.*DeltaPhi)'*G)';
dVda=ep.*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(TeleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
2 The table of input parameters and the code of Example 2
|
nelx |
nely |
volfrac |
ep |
dsp |
kmin |
|
60 |
60 |
0.5 |
0.01 |
3 |
0.035 |
function TOCSPRBF_MMA(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=60;volfrac=0.5;ep=0.01;nele = nelx*nely;
%% LEVEL SET FUNCTION INITIALIZATION水平集函数初始化
r=nely*0.05;
hX = nelx*[repmat([1/6,1/3,1/2,2/3,5/6],1,5)];
hY = nely*[kron([1/6,1/3,1/2,2/3,5/6],ones(1,5))];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION径向基函数初始化
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=10000;d=0;m=1;%MMA参数
%% THERMAL FINITE ELEMENT ANALYSIS PREPARATION
TA1 = [ 8 -2; -2 8];
TA2 = [-4 -2; -2 -4];
KEth = 1/12*[TA1 TA2; TA2 TA1];
TeleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
TeleNode = repmat(TeleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
TedofMat = TeleNode;
TiK = reshape(kron(TedofMat,ones(4,1))',16*nele,1);
TjK = reshape(kron(TedofMat,ones(1,4))',16*nele,1);
Tmaxdof = (nely+1)*(nelx+1);
% DEFINE LOADS AND SUPPORTS (HEATED BOTTOM EDGE)
k0 = 1; kmin = 0.03;
T0 = 0;
Qnode =[(nely+1)*nelx/3+nely/2+1,(nely+1)*nelx/2+nely/3+1,(nely+1)*nelx/2+nely*2/3+1,(nely+1)*nelx*2/3+nely/2+1];
Q = sparse(Qnode,1,10,nNode,1);
T = ones(Tmaxdof,1)*T0;
Tcdofs = [1 nely+1 (nely+1)*nelx+1 (nely+1)*(nelx+1)];
Tfdofs = setdiff(1:Tmaxdof,Tcdofs);
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 100; nRelax = 10;
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(TeleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(TeleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(TeleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(TeleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
TsK = reshape(KEth(:)*(kmin+eleVol'*(k0-kmin)),16*nelx*nely,1);
Kth = sparse(TiK,TjK,TsK); Kth = (Kth+Kth')/2;
T(Tfdofs) = Kth(Tfdofs,Tfdofs)\(Q(Tfdofs)-Kth(Tfdofs,Tcdofs)*T(Tcdofs));
TeleComp = 1/2*sum((T(TedofMat)*KEth).*T(TedofMat),2).*(kmin+eleVol*(k0-kmin));
comp(iT) = sum(TeleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
TeleComp = reshape(TeleComp,nely,nelx);
TeleCompLR = [TeleComp(:,1),TeleComp]+[TeleComp,TeleComp(:,end)];
TnodeComp = ([TeleCompLR;TeleCompLR(end,:)]+[TeleCompLR(1,:);TeleCompLR])/4;
%% MMA准备
DeltaPhi=reshape(DeltaPhi,nNode,1);
TnodeComp=reshape(TnodeComp,nNode,1);
dJda=-((TnodeComp.*DeltaPhi)'*G)';
dVda=ep.*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(TeleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
3 The table of input parameters and the code of Example 3
|
nelx |
nely |
volfrac |
ep |
dsp |
kmin |
|
60 |
60 |
0.5 |
0.01 |
3 |
0.05 |
function TOCSPRBF_MMA(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=60;volfrac=0.5;ep=0.01;nele = nelx*nely;
%% LEVEL SET FUNCTION INITIALIZATION水平集函数初始化
r=nely*0.05;
hX = nelx*[repmat([1/7,2/7,3/7,4/7,5/7,6/7],1,6)];
hY = nely*[kron([1/7,2/7,3/7,4/7,5/7,6/7],ones(1,6))];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION径向基函数初始化
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=10000;d=0;m=1;%MMA参数
%% THERMAL FINITE ELEMENT ANALYSIS PREPARATION
TA1 = [ 8 -2; -2 8];
TA2 = [-4 -2; -2 -4];
KEth = 1/12*[TA1 TA2; TA2 TA1];
TeleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
TeleNode = repmat(TeleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
TedofMat = TeleNode;
TiK = reshape(kron(TedofMat,ones(4,1))',16*nele,1);
TjK = reshape(kron(TedofMat,ones(1,4))',16*nele,1);
Tmaxdof = (nely+1)*(nelx+1);
% DEFINE LOADS AND SUPPORTS (HEATED BOTTOM EDGE)
k0 = 1; kmin = 0.05;
T0 = 0;
Q = sparse([],[],1,Tmaxdof,1);
T = ones(Tmaxdof,1)*T0;
TBC1 = 1:nely+1:Tmaxdof;
TBC2 = nely+1:nely+1:Tmaxdof;
TBC3 = 1:nely+1;
TBC4 = (nely+1)*nelx+1:(nely+1)*(nelx+1);
Tcdofs = [TBC1 TBC2 TBC3 TBC4];
Tfdofs = setdiff(1:Tmaxdof,Tcdofs);
Q(Tfdofs)=1;
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 100; nRelax = 10;
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(TeleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(TeleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(TeleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(TeleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
TsK = reshape(KEth(:)*(kmin+eleVol'*(k0-kmin)),16*nelx*nely,1);
Kth = sparse(TiK,TjK,TsK); Kth = (Kth+Kth')/2;
T(Tfdofs) = Kth(Tfdofs,Tfdofs)\(Q(Tfdofs)-Kth(Tfdofs,Tcdofs)*T(Tcdofs));
TeleComp = 1/2*sum((T(TedofMat)*KEth).*T(TedofMat),2).*(kmin+eleVol*(k0-kmin));
comp(iT) = sum(TeleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
TeleComp = reshape(TeleComp,nely,nelx);
TeleCompLR = [TeleComp(:,1),TeleComp]+[TeleComp,TeleComp(:,end)];
TnodeComp = ([TeleCompLR;TeleCompLR(end,:)]+[TeleCompLR(1,:);TeleCompLR])/4;
%% MMA准备
DeltaPhi=reshape(DeltaPhi,nNode,1);
TnodeComp=reshape(TnodeComp,nNode,1);
dJda=-((TnodeComp.*DeltaPhi)'*G)';
dVda=ep.*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(TeleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
4 The table of input parameters and the code of Example 4
|
nelx |
nely |
volfrac |
ep |
dsp |
kmin |
|
60 |
60 |
0.5 |
0.02 |
3 |
0.06 |
function TOCSPRBF_MMA(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=60;volfrac=0.5;ep=0.01;nele = nelx*nely;
%% LEVEL SET FUNCTION INITIALIZATION水平集函数初始化
r=nely*0.05;
hX = nelx*[repmat([1/6,1/3,1/2,2/3,5/6],1,5)];
hY = nely*[kron([1/6,1/3,1/2,2/3,5/6],ones(1,5))];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION径向基函数初始化
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=10000;d=0;m=1;%MMA参数
%% THERMAL FINITE ELEMENT ANALYSIS PREPARATION
TA1 = [ 8 -2; -2 8];
TA2 = [-4 -2; -2 -4];
KEth = 1/12*[TA1 TA2; TA2 TA1];
TeleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
TeleNode = repmat(TeleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
TedofMat = TeleNode;
TiK = reshape(kron(TedofMat,ones(4,1))',16*nele,1);
TjK = reshape(kron(TedofMat,ones(1,4))',16*nele,1);
Tmaxdof = (nely+1)*(nelx+1);
% DEFINE LOADS AND SUPPORTS (HEATED BOTTOM EDGE)
k0 = 1; kmin = 0.06;
T0 = 0;
Q = sparse([],[],0,Tmaxdof,1);
T = ones(Tmaxdof,1)*T0;
Tcdofs = [1 nely+1 (nely+1)*nelx+1 (nely+1)*(nelx+1)];
Tfdofs = setdiff(1:Tmaxdof,Tcdofs);
Q(Tfdofs)=1;
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 100; nRelax = 10;
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(TeleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(TeleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(TeleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(TeleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
TsK = reshape(KEth(:)*(kmin+eleVol'*(k0-kmin)),16*nelx*nely,1);
Kth = sparse(TiK,TjK,TsK); Kth = (Kth+Kth')/2;
T(Tfdofs) = Kth(Tfdofs,Tfdofs)\(Q(Tfdofs)-Kth(Tfdofs,Tcdofs)*T(Tcdofs));
TeleComp = 1/2*sum((T(TedofMat)*KEth).*T(TedofMat),2).*(kmin+eleVol*(k0-kmin));
comp(iT) = sum(TeleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
TeleComp = reshape(TeleComp,nely,nelx);
TeleCompLR = [TeleComp(:,1),TeleComp]+[TeleComp,TeleComp(:,end)];
TnodeComp = ([TeleCompLR;TeleCompLR(end,:)]+[TeleCompLR(1,:);TeleCompLR])/4;
%% MMA准备
DeltaPhi=reshape(DeltaPhi,nNode,1);
TnodeComp=reshape(TnodeComp,nNode,1);
dJda=-((TnodeComp.*DeltaPhi)'*G)';
dVda=ep.*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(TeleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
5 The table of input parameters and the code of Example 5
|
nelx |
nely |
volfrac |
ep |
dsp |
kmin |
|
60 |
60 |
0.5 |
0.08 |
3 |
0.035 |
function TOCSPRBF_MMA(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=60;volfrac=0.5;ep=0.01;nele = nelx*nely;
%% LEVEL SET FUNCTION INITIALIZATION水平集函数初始化
r=nely*0.05;
hX = nelx*[repmat([1/7,2/7,3/7,4/7,5/7,6/7],1,6)];
hY = nely*[kron([1/7,2/7,3/7,4/7,5/7,6/7],ones(1,6))];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION径向基函数初始化
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=10000;d=0;m=1;%MMA参数
%% THERMAL FINITE ELEMENT ANALYSIS PREPARATION
TA1 = [ 8 -2; -2 8];
TA2 = [-4 -2; -2 -4];
KEth = 1/12*[TA1 TA2; TA2 TA1];
TeleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
TeleNode = repmat(TeleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
TedofMat = TeleNode;
TiK = reshape(kron(TedofMat,ones(4,1))',16*nele,1);
TjK = reshape(kron(TedofMat,ones(1,4))',16*nele,1);
Tmaxdof = (nely+1)*(nelx+1);
% DEFINE LOADS AND SUPPORTS (HEATED BOTTOM EDGE)
k0 = 1; kmin = 0.035;
T0 = 0;
Q = sparse([],[],0,Tmaxdof,1);
T = ones(Tmaxdof,1)*T0;
Tcdofs = [nely/2-1,nely/2,nely/2+1,nely/2+2,nely/2+3];
Tfdofs = setdiff(1:Tmaxdof,Tcdofs);
Q(Tfdofs)=1;
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 100; nRelax = 10;
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(TeleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(TeleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(TeleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(TeleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
TsK = reshape(KEth(:)*(kmin+eleVol'*(k0-kmin)),16*nelx*nely,1);
Kth = sparse(TiK,TjK,TsK); Kth = (Kth+Kth')/2;
T(Tfdofs) = Kth(Tfdofs,Tfdofs)\(Q(Tfdofs)-Kth(Tfdofs,Tcdofs)*T(Tcdofs));
TeleComp = 1/2*sum((T(TedofMat)*KEth).*T(TedofMat),2).*(kmin+eleVol*(k0-kmin));
comp(iT) = sum(TeleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
TeleComp = reshape(TeleComp,nely,nelx);
TeleCompLR = [TeleComp(:,1),TeleComp]+[TeleComp,TeleComp(:,end)];
TnodeComp = ([TeleCompLR;TeleCompLR(end,:)]+[TeleCompLR(1,:);TeleCompLR])/4;
%% MMA准备
DeltaPhi=reshape(DeltaPhi,nNode,1);
TnodeComp=reshape(TnodeComp,nNode,1);
dJda=-((TnodeComp.*DeltaPhi)'*G)';
dVda=ep.*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(TeleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
（The above contents were uploaded on 22 Jan 2025）
***************************************************************************************************************************************************************************************
***************************************************************************************************************************************************************************************
***************************************************************************************************************************************************************************************
（The following content were uploaded on 11 Jan 2025）
The codes and data of the article "A parameterized level set method for structural topology optimization using the approximate re-initialization scheme"
1 The table of input parameters and the code of Example 1
|
nelx |
nely |
volfrac |
ep |
dsp |
Alphamin |
Alphamax |
|
60 |
30 |
0.5 |
0.02 |
3 |
-3 |
3 |
function TOCSPRBF_MMA1(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');%MMA file location
tic;
nelx=60;nely=30;volfrac=0.5;ep=0.02;
%% LEVEL SET FUNCTION INITIALIZATION
r=nely*0.1;%RADIUS OF INITIAL HOLES
hX = nelx*[repmat([1/6,5/6],1,3),repmat([0,1/3,2/3,1],1,2),1/2];
hY = nely*[kron([0,1/2,1],ones(1,2)),kron([1/4,3/4],ones(1,4)),1/2];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin =-3.*ones(nNode,1);
Alphamax=3.*ones(nNode,1);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=100000;d=0;m=1;%MMA
%% FINITE ELEMENT ANALYSIS PREPARATION
E0 = 1; Emin = 1e-9; nu = 0.3; %MATERIAL PROPERTIESE0
A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
eleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
eleNode = repmat(eleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
edofMat = kron(eleNode,[2,2])+repmat([-1,0],nelx*nely,4);
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
%% BOUNDARY CONDITION DEFINITION
F = sparse(2*((nely+1)*nelx+ceil(nely/2)+1),1,-1,2*nNode,1);
fixeddofs = 1:1:2*(nely+1);
freedofs = setdiff(1:2*nNode,fixeddofs);
U = zeros(2*nNode,1);
%% ITERATION OPTIMIZATION
nLoop = 200; nRelax=10;%nLoop
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(eleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
sK = reshape(KE(:)*(Emin+eleVol'*(E0-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs,1) = K(freedofs,freedofs)\F(freedofs,1);
eleComp = sum((U(edofMat)*KE).*U(edofMat),2).*(Emin+eleVol*(E0-Emin));
comp(iT) = sum(eleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
eleComp = reshape(eleComp,nely,nelx);
eleCompLR = [eleComp(:,1),eleComp]+[eleComp,eleComp(:,end)];
nodeComp = ([eleCompLR;eleCompLR(end,:)]+[eleCompLR(1,:);eleCompLR])/4;
%% MMA
DeltaPhi=reshape(DeltaPhi,nNode,1);
nodeComp=reshape(nodeComp,nNode,1);
dJda=-((nodeComp.*DeltaPhi)'*G)';
dVda=ep*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(eleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
2 The table of input parameters and the code of Example 2
|
nelx |
nely |
volfrac |
ep |
dsp |
Alphamin |
Alphamax |
|
120 |
30 |
0.5 |
0.02 |
3 |
-3 |
3 |
function TOCSPRBF_MMA2(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');%MMA file location
tic;
nelx=120;nely=30;volfrac=0.5;ep=0.02;
%% LEVEL SET FUNCTION INITIALIZATION
r=nely*0.05;%RADIUS OF INITIAL HOLES
hX = nelx*[repmat([1/6,5/6],1,3),repmat([0,1/3,2/3,1],1,2),1/2];
hY = nely*[kron([0,1/2,1],ones(1,2)),kron([1/4,3/4],ones(1,4)),1/2];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=-3.*ones(nNode,1);
Alphamax=3.*ones(nNode,1);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=100000;d=0;m=1;%MMA
%% FINITE ELEMENT ANALYSIS PREPARATION
E0 = 1; Emin = 1e-9; nu = 0.3; %MATERIAL PROPERTIESE0
A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
eleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
eleNode = repmat(eleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
edofMat = kron(eleNode,[2,2])+repmat([-1,0],nelx*nely,4);
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
%% BOUNDARY CONDITION DEFINITION 边界条件定义
F = sparse(2*((nely+1)*(ceil(nelx/2)+1)),1,-1,2*nNode,1);%NODAL LOADS
fixeddofs =[1,2,2*(nely+1)*nelx+2];%DISPLACEMENT CONSTRAINTS%桥
freedofs = setdiff(1:2*nNode,fixeddofs);
U = zeros(2*nNode,1);
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 200; nRelax=10;%nLoop 优化的最大迭代次数
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(eleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
sK = reshape(KE(:)*(Emin+eleVol'*(E0-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs,1) = K(freedofs,freedofs)\F(freedofs,1);
eleComp = sum((U(edofMat)*KE).*U(edofMat),2).*(Emin+eleVol*(E0-Emin));
comp(iT) = sum(eleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
eleComp = reshape(eleComp,nely,nelx);
eleCompLR = [eleComp(:,1),eleComp]+[eleComp,eleComp(:,end)];
nodeComp = ([eleCompLR;eleCompLR(end,:)]+[eleCompLR(1,:);eleCompLR])/4;
%% MMA
DeltaPhi=reshape(DeltaPhi,nNode,1);
nodeComp=reshape(nodeComp,nNode,1);
dJda=-((nodeComp.*DeltaPhi)'*G)';
dVda=ep*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(eleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
3 The table of input parameters and the code of Example 3
|
nelx |
nely |
volfrac |
ep |
dsp |
Alphamin |
Alphamax |
|
30 |
60 |
0.3 |
0.01 |
3 |
-3 |
3 |
function TOCSPRBF_MMA3(nelx,nely,volfrac,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=30;nely=60;volfrac=0.3;ep=0.01;
%% LEVEL SET FUNCTION INITIALIZATION
r=nely*0.05;%RADIUS OF INITIAL HOLES
hX = nelx*[repmat([1/6,5/6],1,3),repmat([0,1/3,2/3,1],1,2),1/2];
hY = nely*[kron([0,1/2,1],ones(1,2)),kron([1/4,3/4],ones(1,4)),1/2];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi = max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r,[],3)));
%% RADIAL BASIS FUNCTION INITIALIZATION
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha = G\[Phi(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=-3.*ones(nNode,1);
Alphamax=3.*ones(nNode,1);
low=Alphamin;
upp=Alphamax;
a0=1;a=0;c=100000;d=0;m=1;%MMA
%% FINITE ELEMENT ANALYSIS PREPARATION
E0 = 1; Emin = 1e-9; nu = 0.3; %MATERIAL PROPERTIESE0
A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
eleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
eleNode = repmat(eleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
edofMat = kron(eleNode,[2,2])+repmat([-1,0],nelx*nely,4);
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
%% BOUNDARY CONDITION DEFINITION 边界条件定义
F = sparse(2*((nely+1)*nelx+ceil(nely/2)+1),1,-1,2*nNode,1);%NODAL LOADS
fixeddofs = 1:1:2*(nely+1);%DISPLACEMENT CONSTRAINTS%桥
freedofs = setdiff(1:2*nNode,fixeddofs);
U = zeros(2*nNode,1);
%% ITERATION OPTIMIZATION 主优化循环
nLoop = 200; nRelax=10;%nLoop
delta = 10;
comp = zeros(nLoop,1); vol = zeros(nLoop,1); change=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi = (1-s(:)).*(1-t(:))/4*Phi(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi(eleNode(:,4))';
eleVol = sum(tmpPhi>=0,1)'/numel(s);
vol(iT) = sum(eleVol)/(nelx*nely);
sK = reshape(KE(:)*(Emin+eleVol'*(E0-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs,1) = K(freedofs,freedofs)\F(freedofs,1);
eleComp = sum((U(edofMat)*KE).*U(edofMat),2).*(Emin+eleVol*(E0-Emin));
comp(iT) = sum(eleComp);
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
fval=vol(iT)-volfrac;
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),vol(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,[0,0]);
colormap([0,0,0]); set(gcf,'color','w'); axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 && ...
all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta = (abs(Phi(:))<=delta);
DeltaPhi = zeros(size(Phi));
DeltaPhi(indexDelta) = 0.75/delta*(1-Phi(indexDelta).^2/delta^2);
eleComp = reshape(eleComp,nely,nelx);
eleCompLR = [eleComp(:,1),eleComp]+[eleComp,eleComp(:,end)];
nodeComp = ([eleCompLR;eleCompLR(end,:)]+[eleCompLR(1,:);eleCompLR])/4;
%% MMA
DeltaPhi=reshape(DeltaPhi,nNode,1);
nodeComp=reshape(nodeComp,nNode,1);
dJda=-((nodeComp.*DeltaPhi)'*G)';
dVda=ep*(DeltaPhi'*G);
[Alphamma,~,~,~,~,~,~,~,~,low,upp] = ...
mmasub(m,nNode,iT,Alpha,Alphamin,Alphamax,Alphaold1,Alphaold2,comp(iT),dJda,fval,dVda,low,upp,a0,a,c,d);
Alpha = Alphamma/mean(gradPhi(unique(eleNode((eleVol<1 & eleVol>0),:))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi = reshape(G*Alpha,nely+1,nelx+1);
end
end
4 The table of input parameters and the code of Example 4
|
nelx |
nely |
Vmax2 |
Vmax3 |
ep |
dsp |
Alphamin |
Alphamax |
|
60 |
30 |
0.5 |
0.3 |
0.01 |
3 |
-3 |
3 |
function OCSPRBF_MMA4(nelx,nely,,Vmax2,Vmax3,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=30;ep=0.01;
Vmax2=0.5;
Vmax3=0.3;
%% LEVEL SET FUNCTION INITIALIZATION
r1 = nely*0.05;
r2=nely*0.13;%RADIUS OF INITIAL HOLES
hX = nelx*[repmat([1/6,5/6],1,3),repmat([0,1/3,2/3,1],1,2),1/2];
hY = nely*[kron([0,1/2,1],ones(1,2)),kron([1/4,3/4],ones(1,4)),1/2];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi2= max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r1,[],3)));
Phi3 =-max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r2,[],3)));
Phi =zeros(nely+1,nelx+1);
index3=Phi2>0;
Phi(index3)=1;
index4=(Phi2>0 & Phi3<0);
Phi(index4)=2;
%% RADIAL BASIS FUNCTION INITIALIZATION
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha(:,1) = G\[Phi2(:)];
Alpha(:,2) = G\[Phi3(:)];
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=-2.9*ones(nNode,2);
Alphamax=2.9*ones(nNode,2);
low=Alphamin;
upp=Alphamax;
a0=1;a=[0,0]';c=[100000,100000]';d=[0,0]';m=1;%MMA
%% FINITE ELEMENT ANALYSIS PREPARATION
E1 = 1; E2=3; Emin = 1e-9; nu = 0.3; %MATERIAL PROPERTIESE0
A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
eleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
eleNode = repmat(eleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
edofMat = kron(eleNode,[2,2])+repmat([-1,0],nelx*nely,4);
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
%% BOUNDARY CONDITION DEFINITION
F = sparse(2*((nely+1)*nelx+ceil(nely/2)+1),1,-1,2*nNode,1);
fixeddofs = 1:1:2*(nely+1);
freedofs = setdiff(1:2*nNode,fixeddofs);
U = zeros(2*nNode,1);
%% ITERATION OPTIMIZATION
nLoop = 200; nRelax = 30;%nLoop
delta = 10;
comp = zeros(nLoop,1); vol2=zeros(nLoop,1);vol3=zeros(nLoop,1);
change=zeros(nLoop,1);f=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi2 = (1-s(:)).*(1-t(:))/4*Phi2(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi2(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi2(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi2(eleNode(:,4))';
eleVol2 = sum(tmpPhi2>=0,1)'/numel(s);
tmpPhi3 = (1-s(:)).*(1-t(:))/4*Phi3(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi3(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi3(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi3(eleNode(:,4))';
eleVol3 = sum(tmpPhi3<=0,1)'/numel(s).*eleVol2;
H2=reshape(derac(Phi2,delta),(nely+1)*(nelx+1),1);
H3=reshape(derac(Phi3,delta),(nely+1)*(nelx+1),1);
vol2(iT) = sum(eleVol2)/(nelx*nely);
vol3(iT) = sum(eleVol3)/(nelx*nely);
V2(iT)=vol2(iT)-vol3(iT);
% FINITE ELEMENT ANALYSIS
sK = reshape(KE(:)*(Emin+(eleVol2-eleVol3)'*(E1-Emin)+eleVol3'*(E2-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs,1) = K(freedofs,freedofs)\F(freedofs,1);
SED=sum((U(edofMat)*KE).*U(edofMat),2);
eleComp = SED.*(Emin+eleVol2*(E1-Emin)+eleVol3*(E2-E1));
eleComp2=SED.*(Emin+eleVol2*(E1-Emin));
eleComp3 = SED.*(Emin+eleVol3*(E1-E2));
comp(iT) = sum(eleComp);
fval2=(vol2(iT)-Vmax2);
fval3=(vol3(iT)-Vmax3);
fval=[fval2;fval3];
if iT==1
change(iT)=0;
else
change(iT)=comp(iT)-comp(iT-1);
end
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol1:%f, Vol2:%f, change:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),V2(iT),vol3(iT),change(iT),toc-time0,toc]);
figure(1); contourf(Phi,2);
colormap([1 1 1;1 0 0;0 1 0]); %set(gcf,'color','w');
axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot([1:iT],vol2(1:iT),'-',[1:iT],vol3(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax &&abs(vol2(iT)-Vmax2)/Vmax2<1e-3 ...
&& abs(vol3(iT)-Vmax3)/Vmax3<1e-3...
&&all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta2 = (abs(Phi2(:))<=delta);
DeltaPhi2 = zeros(size(Phi2));
DeltaPhi2(indexDelta2) = 0.75/delta*(1-Phi2(indexDelta2).^2/delta^2);
indexDelta3 = (abs(Phi3(:))<=delta);
DeltaPhi3 = zeros(size(Phi3));
DeltaPhi3(indexDelta3) = 0.75/delta*(1-Phi3(indexDelta3).^2/delta^2);
eleComp2 = reshape(eleComp,nely,nelx);
eleCompLR2 = [eleComp2(:,1),eleComp2]+[eleComp2,eleComp2(:,end)];
nodeComp2 = ([eleCompLR2;eleCompLR2(end,:)]+[eleCompLR2(1,:);eleCompLR2])/4;
eleComp3 = reshape(eleComp3,nely,nelx);
eleCompLR3 = [eleComp3(:,1),eleComp3]+[eleComp3,eleComp3(:,end)];
nodeComp3 = ([eleCompLR3;eleCompLR3(end,:)]+[eleCompLR3(1,:);eleCompLR3])/4;
%% MMA
DeltaPhi2=reshape(DeltaPhi2,nNode,1);
DeltaPhi3=reshape(DeltaPhi3,nNode,1);
nodeComp2=reshape(nodeComp2,nNode,1);
nodeComp3=reshape(nodeComp3,nNode,1);
dJda(:,1)=-((nodeComp2.*DeltaPhi2)'*G);
dJda(:,2)=-((nodeComp3.*DeltaPhi3)'*G);
dVda(:,1)=((DeltaPhi2)'*G);
dVda(:,2)=-((H2.*DeltaPhi3)'*G);
dVda=ep.*dVda;
[Alphamma(:,1),~,~,~,~,~,~,~,~,low(:,1),upp(:,1)] = ...
mmasub(1,nNode,iT,Alpha(:,1),Alphamin(:,1),Alphamax(:,1),Alphaold1(:,1),Alphaold2(:,1),comp(iT),dJda(:,1),fval2,dVda(:,1)',low(:,1),upp(:,1),a0,a(1),c(1),d(1));
[Alphamma(:,2),~,~,~,~,~,~,~,~,low(:,2),upp(:,2)] = ...
mmasub(1,nNode,iT,Alpha(:,2),Alphamin(:,2),Alphamax(:,2),Alphaold1(:,2),Alphaold2(:,2),comp(iT),dJda(:,2),fval3,dVda(:,2)',low(:,2),upp(:,2),a0,a(2),c(2),d(2));
Alpha(:,1) = Alphamma(:,1)/mean(gradPhi(unique(eleNode((eleVol2<1 & eleVol2>0),1))));
Alpha(:,2) = Alphamma(:,2)/mean(gradPhi(unique(eleNode((eleVol3<1 & eleVol3>0),2))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi2 = G*Alpha(:,1);
Phi3 = G*Alpha(:,2);
Phi=zeros(nNode,1);
index3=(Phi2(:)>0);
Phi(index3)=1;
index4=(Phi2(:)>0 & Phi3(:)<0);
Phi(index4)=2;
Phi=reshape(Phi,nely+1,nelx+1);
end
end
function [heaviside]=derac(phi,dphi)
% Dirac delta & heaviside functions
m=numel(phi);
delta=0.0001.*ones(m,1);
heaviside=ones(m,1);
for i=1:m
if abs(phi(i))<=dphi
delta(i)=3*(1-0.0001)/(4*dphi)*(1-(abs(phi(i)))^2/(dphi^2));
heaviside(i)=3/4*(phi(i)/dphi-(phi(i)^3)/(3*dphi^3))+0.5;
elseif phi(i)<=-dphi
heaviside(i)=0.00005;
end
end
end
5 The table of input parameters and the code of Example 5
|
nelx |
nely |
Vmax1 |
Vmax2 |
Vmax13 |
ep |
dsp |
Alphamin |
Alphamax |
|
30 |
60 |
0.5 |
0.4 |
0.3 |
0.02 |
3 |
-3 |
3 |
function TOCSPRBF_MMA5(nelx,nely,Vmax1,Vmax2,Vmax3,ep)
addpath('C:\Users180\Documents\MATLAB\GCMMA-MMA-code-1.5\');tic;
nelx=60;nely=30;ep=0.02;Vmax1=0.5;Vmax2=0.4;Vmax3=0.3;
%% LEVEL SET FUNCTION INITIALIZATION
r1 =nely*0.05;
r2=nely*0.1;
r3=nely*0.15;%RADIUS OF INITIAL HOLES
hX = nelx*[repmat([1/6,5/6],1,3),repmat([0,1/3,2/3,1],1,2),1/2];
hY = nely*[kron([0,1/2,1],ones(1,2)),kron([1/4,3/4],ones(1,4)),1/2];
[X,Y] = meshgrid(0:1:nelx,0:1:nely);
dX = bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
dY = bsxfun(@minus,repmat(Y,[1,1,numel(hY)]),reshape(hY,1,1,numel(hY)));
Phi1= max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r1,[],3)));
Phi2= -max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r2,[],3)));
Phi3= -max(-3,min(3,min(sqrt(dX.^2+dY.^2)-r3,[],3)));
Phi =zeros(nely+1,nelx+1);
index1=Phi1>0;
Phi(index1)=1;
index2=(Phi1>0 & Phi2<0);%..
Phi(index2)=2;
index3=(Phi1>0 & Phi2<0 & Phi3<0);
Phi(index3)=3;
%% RADIAL BASIS FUNCTION INITIALIZATION
cRBF = 1e-4;%RBF PARAMETER
nNode = (nely+1)*(nelx+1);
Ax = bsxfun(@minus,X(:),X(:)');
Ay = bsxfun(@minus,Y(:),Y(:)');
dsp=3;
rCS = sqrt(Ax.^2+Ay.^2+cRBF^2)/dsp;%CSRBF
G=(max(0,1-rCS).^4).*(4*rCS+1);
pGpX = (max(0,1-rCS).^3).*(-20*rCS).*(Ax./rCS/dsp^2);
pGpY = (max(0,1-rCS).^3).*(-20*rCS).*(Ay./rCS/dsp^2);
Alpha(:,1) = G\[Phi1(:)];
Alpha(:,2) = G\[Phi2(:)];
Alpha(:,3) = G\[Phi3(:)];
% mma
Alphaold1=Alpha;
Alphaold2=Alpha;
Alphamin=-3*ones(nNode,3);
Alphamax=3*ones(nNode,3);
low=Alphamin;
upp=Alphamax;
a0=1;a=[0,0,0]';c=[100000,100000,100000]';d=[0,0,0]';m=1;%MMA
%% FINITE ELEMENT ANALYSIS PREPARATION
E1 = 1; E2=3; E3 = 9; Emin = 1e-9; nu = 0.3; %MATERIAL PROPERTIESE0
A11 = [12 3 -6 -3; 3 12 3 0; -6 3 12 -3; -3 0 -3 12];
A12 = [-6 -3 0 3; -3 -6 -3 -6; 0 -3 -6 3; 3 -6 3 -6];
B11 = [-4 3 -2 9; 3 -4 -9 4; -2 -9 -4 -3; 9 4 -3 -4];
B12 = [ 2 -3 4 -9; -3 2 9 -2; 4 9 2 3; -9 -2 3 2];
KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
eleN1 = repmat((1:nely)',1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
eleNode = repmat(eleN1(:),1,4)+repmat([0,nely+[1,2],1],nelx*nely,1);
edofMat = kron(eleNode,[2,2])+repmat([-1,0],nelx*nely,4);
iK = reshape(kron(edofMat,ones(8,1))',64*nelx*nely,1);
jK = reshape(kron(edofMat,ones(1,8))',64*nelx*nely,1);
%% BOUNDARY CONDITION DEFINITION
F = sparse(2*((nely+1)*nelx+ceil(nely/2)+1),1,-1,2*nNode,1);
fixeddofs = 1:1:2*(nely+1);
freedofs = setdiff(1:2*nNode,fixeddofs);
U = zeros(2*nNode,1);
%% ITERATION OPTIMIZATION
nLoop =100; nRelax = 30;%nLoop
delta = 10;
comp = zeros(nLoop,1); vol1 = zeros(nLoop,1);vol2=zeros(nLoop,1);vol3=zeros(nLoop,1);
change=zeros(nLoop,1);V1=zeros(nLoop,1);V2=zeros(nLoop,1);
for iT = 1:nLoop
time0 = toc;
%% FINITE ELEMENT ANALYSIS
[s,t] = meshgrid(-1:0.1:1,-1:0.1:1);
tmpPhi1 = (1-s(:)).*(1-t(:))/4*Phi1(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi1(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi1(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi1(eleNode(:,4))';
eleVol1 = sum(tmpPhi1>=0,1)'/numel(s);
tmpPhi2 = (1-s(:)).*(1-t(:))/4*Phi2(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi2(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi2(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi2(eleNode(:,4))';
eleVol2 = sum(tmpPhi2<=0,1)'/numel(s).*eleVol1;
tmpPhi3 = (1-s(:)).*(1-t(:))/4*Phi3(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*...
Phi3(eleNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi3(eleNode(:,3))'+...
(1-s(:)).*(1+t(:))/4*Phi3(eleNode(:,4))';
eleVol3 = sum(tmpPhi3<=0,1)'/numel(s).*eleVol2;
H1=reshape(derac(Phi1,delta),(nely+1)*(nelx+1),1);
H2=reshape(derac(Phi2,delta),(nely+1)*(nelx+1),1);
H3=reshape(derac(Phi3,delta),(nely+1)*(nelx+1),1);
vol1(iT) = sum(eleVol1)/(nelx*nely);
vol2(iT) = sum(eleVol2)/(nelx*nely);
vol3(iT) = sum(eleVol3)/(nelx*nely);
sK = reshape(KE(:)*(Emin+(eleVol1-eleVol2)'*(E1-Emin)+(eleVol2-eleVol3)'*(E2-Emin)+eleVol3'*(E3-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
U(freedofs,1) = K(freedofs,freedofs)\F(freedofs,1);
SED=sum((U(edofMat)*KE).*U(edofMat),2);
eleComp = SED.*(Emin+eleVol1*(E1-Emin)+eleVol2*(E2-E1)+eleVol3*(E3-E2));
eleComp2=SED.*(Emin+eleVol2*(E1-E2)+eleVol3*(E2-E3));%????
eleComp3=SED.*(Emin+eleVol3*(E2-E3));
comp(iT) = sum(eleComp);
fval1=(vol1(iT)-Vmax1);
fval2=(vol2(iT)-Vmax2);
fval3=(vol3(iT)-Vmax3);
fval=[fval1,fval2,fval3];
V1(iT)=vol1(iT)-vol2(iT);
V2(iT)=vol2(iT)-vol3(iT);
%% DISPLAY RESULTS
fprintf('No.%i, Obj:%f, Vol1:%f, Vol2:%f, Vol3:%f, IterTime:%f, TotalTime:%f\n',[iT,comp(iT),V1(iT),V2(iT),vol3(iT),toc-time0,toc]);
figure(1); contourf(Phi,3);
colormap ([1 1 1;1 0 0; 0 1 0;0 0 1]); %set(gcf,'color','w');
axis equal; axis off;
figure(2); surf(Phi); clim([-12,12]);
axis equal; axis([0,nelx,0,nely,-12,12]); view(3);
figure(3); subplot(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
subplot(2,1,2); plot([1:iT],V1(1:iT),'-',[1:iT],V2(1:iT),'-',[1:iT],vol3(1:iT),'-'); title('Volume fraction');
%% CONVERGENCE CHECK
if iT>nRelax &&abs(vol1(iT)-Vmax1)/Vmax1<1e-3 ...
&& abs(vol2(iT)-Vmax2)/Vmax2<1e-3...
&& abs(vol3(iT)-Vmax3)/Vmax3<1e-3...
&&all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3)
break;
end
%% LEVEL SET FUNCTION EVOLUTION
gradPhi = sqrt((pGpX*Alpha).^2+(pGpY*Alpha).^2);
indexDelta1 = (abs(Phi1(:))<=delta);
DeltaPhi1 = zeros(size(Phi1));
DeltaPhi1(indexDelta1) = 0.75/delta*(1-Phi1(indexDelta1).^2/delta^2);
indexDelta2 = (abs(Phi2(:))<=delta);
DeltaPhi2 = zeros(size(Phi2));
DeltaPhi2(indexDelta2) = 0.75/delta*(1-Phi2(indexDelta2).^2/delta^2);
indexDelta3 = (abs(Phi3(:))<=delta);
DeltaPhi3 = zeros(size(Phi3));
DeltaPhi3(indexDelta3) = 0.75/delta*(1-Phi3(indexDelta3).^2/delta^2);
eleComp1 = reshape(eleComp,nely,nelx);
eleCompLR1 = [eleComp1(:,1),eleComp1]+[eleComp1,eleComp1(:,end)];
nodeComp1 = ([eleCompLR1;eleCompLR1(end,:)]+[eleCompLR1(1,:);eleCompLR1])/4;
eleComp2 = reshape(eleComp2,nely,nelx);
eleCompLR2 = [eleComp2(:,1),eleComp2]+[eleComp2,eleComp2(:,end)];
nodeComp2 = ([eleCompLR2;eleCompLR2(end,:)]+[eleCompLR2(1,:);eleCompLR2])/4;
eleComp3 = reshape(eleComp3,nely,nelx);
eleCompLR3 = [eleComp3(:,1),eleComp3]+[eleComp3,eleComp3(:,end)];
nodeComp3 = ([eleCompLR3;eleCompLR3(end,:)]+[eleCompLR3(1,:);eleCompLR3])/4;
%% MMA
DeltaPhi1=reshape(DeltaPhi1,nNode,1);
DeltaPhi2=reshape(DeltaPhi2,nNode,1);
DeltaPhi3=reshape(DeltaPhi3,nNode,1);
nodeComp1=reshape(nodeComp1,nNode,1);
nodeComp2=reshape(nodeComp2,nNode,1);
nodeComp3=reshape(nodeComp3,nNode,1);
dJda(:,1)=-((nodeComp1.*DeltaPhi1)'*G);
dJda(:,2)=-((nodeComp2.*DeltaPhi1)'*G);
dJda(:,3)=-((nodeComp3.*DeltaPhi3)'*G);
dVda(:,1)=((DeltaPhi1)'*G);
dVda(:,2)=-((H1.*DeltaPhi2)'*G);
dVda(:,3)=-((H1.*(1-H2).*DeltaPhi3)'*G);
dVda=ep.*dVda;
[Alphamma(:,1),~,~,~,~,~,~,~,~,low(:,1),upp(:,1)] = ...
mmasub(1,nNode,iT,Alpha(:,1),Alphamin(:,1),Alphamax(:,1),Alphaold1(:,1),Alphaold2(:,1),comp(iT),dJda(:,1),fval1,dVda(:,1)',low(:,1),upp(:,1),a0,a(1),c(1),d(1));
[Alphamma(:,2),~,~,~,~,~,~,~,~,low(:,2),upp(:,2)] = ...
mmasub(1,nNode,iT,Alpha(:,2),Alphamin(:,2),Alphamax(:,2),Alphaold1(:,2),Alphaold2(:,2),comp(iT),dJda(:,2),fval2,dVda(:,2)',low(:,2),upp(:,2),a0,a(2),c(2),d(2));
[Alphamma(:,3),~,~,~,~,~,~,~,~,low(:,3),upp(:,3)] = ...
mmasub(1,nNode,iT,Alpha(:,3),Alphamin(:,3),Alphamax(:,3),Alphaold1(:,3),Alphaold2(:,3),comp(iT),dJda(:,3),fval3,dVda(:,3)',low(:,3),upp(:,3),a0,a(3),c(3),d(3));
Alpha(:,1) = Alphamma(:,1)/mean(gradPhi(unique(eleNode((eleVol1<1 & eleVol1>0),1))));
Alpha(:,2) = Alphamma(:,2)/mean(gradPhi(unique(eleNode((eleVol1<1 & eleVol1>0),1))));
Alpha(:,3)=Alphamma(:,3)/mean(gradPhi(unique(eleNode((eleVol1<1 & eleVol1>0),1))));
Alphaold2 =Alphaold1;
Alphaold1 = Alpha;
Alphamin=min(Alpha);
Alphamax=max(Alpha);
Phi1 = G*Alpha(:,1);
Phi2 = G*Alpha(:,2);
Phi3 = G*Alpha(:,3);
Phi=zeros(nNode,1);
index4=(Phi1(:)>0);
Phi(index4)=1;
index5=(Phi1(:)>0 & Phi2(:)<0);
Phi(index5)=2;
index6=(Phi1(:)>0 & Phi2(:)<0 & Phi3(:)<0);
Phi(index6)=3;
Phi=reshape(Phi,nely+1,nelx+1);
end
end
function [heaviside]=derac(phi,dphi)
% Dirac delta & heaviside functions
m=numel(phi);
delta=0.0001.*ones(m,1);
heaviside=ones(m,1);
for i=1:m
if abs(phi(i))<=dphi
delta(i)=3*(1-0.0001)/(4*dphi)*(1-(abs(phi(i)))^2/(dphi^2));
heaviside(i)=3/4*(phi(i)/dphi-(phi(i)^3)/(3*dphi^3))+0.5;
elseif phi(i)<=-dphi
heaviside(i)=0.00005;
end
end
end