www.gusucode.com > 基于多边形网格的流动问题拓扑优化设计元代码 > PolyTop.m
%------------------------------ PolyTop ----------------------------------%
% Ref: A Pereira, C Talischi, GH Paulino, IFM Menezes, MS Carvalho %
% "Fluid Flow Topology Optimization in PolyTop: Stability and %
% Computational Implementation", Struct Multidisc Optim, %
% DOI 10.1007/s00158-014-1182-z %
%-------------------------------------------------------------------------%
function [z,V,fem] = PolyTop(fem,opt)
Iter=0; Tol=opt.Tol*(opt.zMax-opt.zMin); Change=2*Tol; z=opt.zIni; P=opt.P;
[E,dEdy,V,dVdy] = opt.MatIntFnc(P*z);
[FigHandle,FigData] = InitialPlot(fem,V);
while (Iter<opt.MaxIter) && (Change>Tol)
Iter = Iter + 1;
%Compute cost functionals and analysis sensitivities
[g,dgdE,dgdV,fem] = ConstraintFnc(fem,E,V,opt.VolFrac);
[f,dfdE,dfdV,fem] = ObjectiveFnc(fem,E,V);
%Compute design sensitivities
dfdz = P'*(dEdy.*dfdE + dVdy.*dfdV);
dgdz = P'*(dEdy.*dgdE + dVdy.*dgdV);
%Update design variable and analysis parameters
[z,Change] = UpdateScheme(dfdz,g,dgdz,z,opt);
[E,dEdy,V,dVdy] = opt.MatIntFnc(P*z);
%Output results
fprintf('It: %i \t Objective: %1.3f\tChange: %1.3f\n',Iter,f,Change);
set(FigHandle,'FaceColor','flat','CData',V(FigData)); drawnow
end
%------------------------------------------------------- OBJECTIVE FUNCTION
function [f,dfdE,dfdV,fem] = ObjectiveFnc(fem,E,V)
[U,F,fem] = FEAnalysis(fem,E);
f = 1/2*dot(F,U);
temp = cumsum(U(fem.iA).*fem.kAalpha.*U(fem.jA));
temp = temp(cumsum(fem.ElemNDofA.^2));
dfdE = 1/2*[temp(1);temp(2:end)-temp(1:end-1)];
dfdV = zeros(size(V));
%------------------------------------------------------ CONSTRAINT FUNCTION
function [g,dgdE,dgdV,fem] = ConstraintFnc(fem,E,V,VolFrac)
if ~isfield(fem,'ElemArea')
fem.ElemArea = zeros(fem.NElem,1);
for el=1:fem.NElem
vx = fem.Node(fem.Element{el},1); vy = fem.Node(fem.Element{el},2);
fem.ElemArea(el) = 0.5*sum(vx.*vy([2:end 1])-vy.*vx([2:end 1]));
end
end
g = sum(fem.ElemArea.*V)/sum(fem.ElemArea)-VolFrac;
dgdE = zeros(size(E));
dgdV = fem.ElemArea/sum(fem.ElemArea);
%----------------------------------------------- OPTIMALITY CRITERIA UPDATE
function [zNew,Change] = UpdateScheme(dfdz,g,dgdz,z0,opt)
zMin=opt.zMin; zMax=opt.zMax;
move=opt.OCMove*(zMax-zMin); eta=opt.OCEta;
l1=0; l2=1e6;
while l2-l1 > 1e-4
lmid = 0.5*(l1+l2);
B = -(dfdz./dgdz)/lmid;
zCnd = zMin+(z0-zMin).*B.^eta;
zNew = max(max(min(min(zCnd,z0+move),zMax),z0-move),zMin);
if (g+dgdz'*(zNew-z0)>0), l1 = lmid;
else l2 = lmid; end
end
Change = max(abs(zNew-z0))/(zMax-zMin);
%-------------------------------------------------------------- FE-ANALYSIS
function [U,F,fem] = FEAnalysis(fem,E)
if ~isfield(fem,'kAmu')
fem.ElemNDofA = 2*cellfun(@length,fem.Element);%% 形成每个单元的自由度索引
fem.iA = zeros(sum(fem.ElemNDofA.^2),1); %% A矩阵的单元组装索引 写成列向量
fem.jA=fem.iA; fem.kAmu=fem.iA; fem.kAalpha=fem.iA; fem.e=fem.iA; %% 类似
fem.iB = zeros(sum(fem.ElemNDofA),1); fem.jB=fem.iB; fem.kB=fem.iB; %% 类似
fem.NDof = 2*fem.NNode+fem.NElem+1; %% 总刚度矩阵中的自由度
indexA = 0; indexB = 0;
if ~isfield(fem,'ShapeFnc'), fem=TabShapeFnc(fem); end
for el = 1:fem.NElem
eNode = fem.Element{el};
eNDof = fem.ElemNDofA(el);
if (el==1 || ~fem.Reg), [Amu_e,Aalpha_e,Be]=LocalK(fem,eNode); end
eDofA = reshape([2*eNode-1;2*eNode],eNDof,1); %% 单元节点对应自由度
I=repmat(eDofA,1,eNDof); J=I'; %% A矩阵的自由度矩阵
fem.iA(indexA+1:indexA+eNDof^2) = I(:); %% 第indexA个单元的行自由度索引
fem.jA(indexA+1:indexA+eNDof^2) = J(:); %% 第indexA个单元的列自由度索引
fem.kAmu(indexA+1:indexA+eNDof^2) = Amu_e(:); %% 第indexA个单元的行列对应Amu_e的值
fem.kAalpha(indexA+1:indexA+eNDof^2) = Aalpha_e(:); %% 用sparse矩阵的命令，可实现自动对应叠加
fem.e(indexA+1:indexA+eNDof^2) = el;
fem.iB(indexB+1:indexB+eNDof) = eDofA(:);
fem.jB(indexB+1:indexB+eNDof) = el;
fem.kB(indexB+1:indexB+eNDof) = Be(:);
indexA = indexA + eNDof^2;
indexB = indexB + eNDof;
end
fem.FixedDofs = zeros(1,size(fem.NodeBC,1));
fem.G = zeros(size(fem.NodeBC,1),1);
for i = 1:size(fem.NodeBC,1)
fem.FixedDofs(i) = 2*(fem.NodeBC(i,1)-1) + fem.NodeBC(i,2);
fem.G(i) = fem.NodeBC(i,3);
end
fem.FreeDofs = setdiff(1:fem.NDof,fem.FixedDofs);
end
A = sparse(fem.iA,fem.jA,fem.kAmu+E(fem.e).*fem.kAalpha);
B = sparse(fem.iB,fem.jB,fem.kB);
Z=zeros(2*fem.NNode,1); O=sparse(fem.NElem,fem.NElem);
K = [A B Z; B' O fem.ElemArea; Z' fem.ElemArea' 0];
K = (K+K')/2;
S = zeros(fem.NDof,1); S(fem.FixedDofs,:)=fem.G;
S(fem.FreeDofs,:) = K(fem.FreeDofs,fem.FreeDofs)\...
(-K(fem.FreeDofs,fem.FixedDofs)*S(fem.FixedDofs,:));
U = S(1:2*fem.NNode);
F = A*U;
%----------------------------------------------- ELEMENT STIFFNESS MATRICES
function [Amu_e,Aalpha_e,Be] = LocalK(fem,eNode)
Cmu = fem.mu0*[2 0 0; 0 2 0; 0 0 1];
nn=length(eNode); Amu_e=zeros(2*nn,2*nn); Aalpha_e=Amu_e; Be=zeros(2*nn,1);
W = fem.ShapeFnc{nn}.W;
for q = 1:length(W) %quadrature loop
dNdxi = fem.ShapeFnc{nn}.dNdxi(:,:,q);
J0 = fem.Node(eNode,:)'*dNdxi;
dNdx = dNdxi/J0; % N=(N1 0 N2 0 N3 0 ... Nnn 0;
% 0 N1 0 N2 0 N3 ... 0 Nnn)
% U=(u1 v1 u2 v2 u3 v3 ... unn vnn)^T
B = zeros(3,2*nn);
B(1,1:2:2*nn) = dNdx(:,1)';
B(2,2:2:2*nn) = dNdx(:,2)';
B(3,1:2:2*nn) = dNdx(:,2)';
B(3,2:2:2*nn) = dNdx(:,1)';
% B=(du/dx 0;0 dv/dy;du/dy+dv/dx)^T
% =(dN1/dx 0 dN2/dx 0 dN3/dx 0 ... dNnn/dx 0
% 0 dN1/dy 0 dN2/dy 0 dN3/dy ... 0 dNnn/dy
% dN1/dy dN1/dx dN2/dy dN2/dx dN3/dy dN3/dx ... dNnn/dy dNnn/dx)
Amu_e = Amu_e + B'*Cmu*B*W(q)*det(J0);
N = fem.ShapeFnc{nn}.N(:,:,q);
Nu = zeros(2,2*nn);
Nu(1,1:2:2*nn) = N(:);
Nu(2,2:2:2*nn) = N(:);
Aalpha_e = Aalpha_e + Nu'*Nu*W(q)*det(J0);
dNudx = reshape(dNdx',1,2*nn);
Be = Be - dNudx'*W(q)*det(J0);
end
%------------------------------------------------- TABULATE SHAPE FUNCTIONS
function fem = TabShapeFnc(fem)
ElemNNode = cellfun(@length,fem.Element); % number of nodes per element
fem.ShapeFnc = cell(max(ElemNNode),1);
for nn = min(ElemNNode):max(ElemNNode)
[W,Q] = PolyQuad(nn);
fem.ShapeFnc{nn}.W = W;
fem.ShapeFnc{nn}.N = zeros(nn,1,size(W,1));
fem.ShapeFnc{nn}.dNdxi = zeros(nn,2,size(W,1));
for q = 1:size(W,1)
[N,dNdxi] = PolyShapeFnc(nn,Q(q,:));
fem.ShapeFnc{nn}.N(:,:,q) = N;
fem.ShapeFnc{nn}.dNdxi(:,:,q) = dNdxi;
end
end
%------------------------------------------------ POLYGONAL SHAPE FUNCTIONS
function [N,dNdxi] = PolyShapeFnc(nn,xi)
N=zeros(nn,1); alpha=zeros(nn,1); dNdxi=zeros(nn,2); dalpha=zeros(nn,2);
sum_alpha=0.0; sum_dalpha=zeros(1,2); A=zeros(nn,1); dA=zeros(nn,2);
[p,Tri] = PolyTrnglt(nn,xi);
for i=1:nn
sctr = Tri(i,:); pT = p(sctr,:);
A(i) = 1/2*det([pT,ones(3,1)]);
dA(i,1) = 1/2*(pT(3,2)-pT(2,2));
dA(i,2) = 1/2*(pT(2,1)-pT(3,1));
end
A=[A(nn,:);A]; dA=[dA(nn,:);dA];
for i=1:nn
alpha(i) = 1/(A(i)*A(i+1));
dalpha(i,1) = -alpha(i)*(dA(i,1)/A(i)+dA(i+1,1)/A(i+1));
dalpha(i,2) = -alpha(i)*(dA(i,2)/A(i)+dA(i+1,2)/A(i+1));
sum_alpha = sum_alpha + alpha(i);
sum_dalpha(1:2) = sum_dalpha(1:2)+dalpha(i,1:2);
end
for i=1:nn
N(i) = alpha(i)/sum_alpha;
dNdxi(i,1:2) = (dalpha(i,1:2)-N(i)*sum_dalpha(1:2))/sum_alpha;
end
%---------------------------------------------------- POLYGON TRIANGULATION
function [p,Tri] = PolyTrnglt(nn,xi)
p = [cos(2*pi*((1:nn))/nn); sin(2*pi*((1:nn))/nn)]';
p = [p; xi];
Tri = zeros(nn,3); Tri(1:nn,1)=nn+1;
Tri(1:nn,2)=1:nn; Tri(1:nn,3)=2:nn+1; Tri(nn,3)=1;
%----------------------------------------------------- POLYGONAL QUADRATURE
function [weight,point] = PolyQuad(nn)
[W,Q]= TriQuad; %integration pnts & wgts for ref. triangle
[p,Tri] = PolyTrnglt(nn,[0 0]); %triangulate from origin
point=zeros(nn*length(W),2); weight=zeros(nn*length(W),1);
for k=1:nn
sctr = Tri(k,:);
for q=1:length(W)
[N,dNds] = TriShapeFnc(Q(q,:)); %compute shape functions
J0 = p(sctr,:)'*dNds;
l = (k-1)*length(W) + q;
point(l,:) = N'*p(sctr,:);
weight(l) = det(J0)*W(q);
end
end
%---------------------------------------------------- TRIANGULAR QUADRATURE
function [weight,point] = TriQuad
point=[1/6,1/6;2/3,1/6;1/6,2/3]; weight=[1/6,1/6,1/6];
%----------------------------------------------- TRIANGULAR SHAPE FUNCTIONS
function [N,dNds] = TriShapeFnc(s)
N=[1-s(1)-s(2);s(1);s(2)]; dNds=[-1,-1;1,0;0,1];
%------------------------------------------------------------- INITIAL PLOT
function [handle,map] = InitialPlot(fem,z0)
Tri = zeros(length([fem.Element{:}])-2*fem.NElem,3);
map = zeros(size(Tri,1),1); index=0;
for el = 1:fem.NElem
for enode = 1:length(fem.Element{el})-2
map(index+1) = el;
Tri(index+1,:) = fem.Element{el}([1,enode+1,enode+2]);
index = index + 1;
end
end
handle = patch('Faces',Tri,'Vertices',fem.Node,'FaceVertexCData',...
z0(map),'FaceColor','flat','EdgeColor','none');
axis equal; axis off; axis tight; colormap(gray); caxis([0 1]);
%-------------------------------------------------------------------------%