[已解决] matlab可以使用.type==1的是什么结构?

[复制链接]
Yintruder 发表于 6 天前
在拓扑优化时,想要跑通一个matlab代码,但是遇到了这样的问题,函数输入有一项‘solver’,但是不知道solver表示什么
  1. function minVpnorm_long(domain,sizex,sizey,helem,penal,rmin,ft,Sytarget,comptarget,...
  2.     betamax,betamaxsimul,solver,filename)
复制代码
在其后,出现了这样的描述:
  1. if solver.type == 4 % Multigrid
复制代码
以及:
  1. solver.cgtol11 = solver.cgtol1;
  2. solver.cgtol22 = solver.cgtol2;
复制代码
这样的数据表示什么意思?
全代码如下:
  1. %%%% 2D TOPOLOGY OPTIMIZATION CODE %%%%
  2. % With stress constraints, p-norm approach
  3. % Single density field at eta=0.5
  4. function minVpnorm_long(domain,sizex,sizey,helem,penal,rmin,ft,Sytarget,comptarget,...
  5.     betamax,betamaxsimul,solver,filename)%betamax=beita_h,betamaxsimul=
  6. tic;
  7. close all;
  8. %% MATERIAL PROPERTIES
  9. Emax = 1;
  10. Emin = Emax*1e-6;
  11. nu = 0.3;
  12. %% PREPARE DOMAIN
  13. switch domain
  14.     case 'biclamped'
  15.         [X,T,i_img,j_img] = generate_biclamped(sizex,sizey,helem,false);
  16.     case 'lbracket'
  17.         [X,T,i_img,j_img] = generate_lbracket(sizex,sizey,helem,false);
  18.     case 'ubracket'
  19.         [X,T,i_img,j_img] = generate_ubracket(sizex,sizey,helem,false);
  20. end
  21. %--------------------------------------------------------------
  22. %% PREPARE FINITE ELEMENT ANALYSIS
  23. A11 = [12  3 -6 -3;  3 12  3  0; -6  3 12 -3; -3  0 -3 12];
  24. A12 = [-6 -3  0  3; -3 -6 -3 -6;  0 -3 -6  3;  3 -6  3 -6];
  25. B11 = [-4  3 -2  9;  3 -4 -9  4; -2 -9 -4 -3;  9  4 -3 -4];
  26. B12 = [ 2 -3  4 -9; -3  2  9 -2;  4  9  2  3; -9 -2  3  2];
  27. KE = 1/(1-nu^2)/24*([A11 A12;A12' A11]+nu*[B11 B12;B12' B11]);
  28. nelem = size(T,1);
  29. nodes = size(X,1);
  30. ndof = 2*nodes;
  31. edofMat = zeros(nelem,8);
  32. edofMat(:,1:2) = [2*T(:,1)-1 2*T(:,1)];
  33. edofMat(:,3:4) = [2*T(:,2)-1 2*T(:,2)];
  34. edofMat(:,5:6) = [2*T(:,3)-1 2*T(:,3)];
  35. edofMat(:,7:8) = [2*T(:,4)-1 2*T(:,4)];
  36. iK = reshape(kron(edofMat,ones(8,1))',64*nelem,1);
  37. jK = reshape(kron(edofMat,ones(1,8))',64*nelem,1);
  38. % Prologation operators
  39. if solver.type == 4 % Multigrid
  40.     Pu = cell(3,1); % Assumes 4 levels
  41.     Xf = X; % Initialize fine grid
  42.     hgrid = helem;
  43.     for i = 1:3
  44.         [Pu{i,1},Xc] = prepcoarse(Xf,hgrid);
  45.         % Replace fine grid with next level
  46.         Xf = Xc;
  47.         hgrid = hgrid*2;
  48.     end
  49. end
  50. U = zeros(ndof,1); LAM = 0*U;
  51. Bmat = 1/helem*[-1/2 0 1/2 0 1/2 0 -1/2 0
  52.     0 -1/2 0 -1/2 0 1/2 0 1/2
  53.     -1/2 -1/2 -1/2 1/2 1/2 1/2 1/2 -1/2];
  54. % The constitutive matrix
  55. Dmat = 1/(1-nu^2)*[ 1 nu 0;nu 1 0;0 0 (1-nu)/2];
  56. %--------------------------------------------------------------
  57. %% DEFINE LOADS AND SUPPORTS
  58. solids = find(T(:,5)==2); % Find solid elements
  59. padding = find(T(:,6)==1); % Find padding elements
  60. solids_pad = intersect(solids,padding); % Find padding solids
  61. loadedelem = setdiff(solids,solids_pad);
  62. nloadedelem = size(loadedelem,1);
  63. switch domain
  64.     case 'ubracket'
  65.         loadeddof = edofMat(loadedelem,1:2:7);
  66.     otherwise
  67.         loadeddof = edofMat(loadedelem,2:2:8);
  68. end
  69. loadeddof = loadeddof(:);
  70. loadmag = -1;
  71. load_per_dof = loadmag/nloadedelem/4;
  72. F = sparse(loadeddof,1,load_per_dof,ndof,1);
  73. % Supports
  74. switch domain
  75.     case 'biclamped'
  76.         % Bi-clamped
  77.         [supnodes,~] = find(X(:,1)==0); % Find nodes with x==0
  78.         supdofs1 = union(2*supnodes-1,2*supnodes);
  79.         [supnodes,~] = find(X(:,1)==sizex); % Find nodes with x==sizex这里也许要删去
  80.         supdofs2 = 2*supnodes-1;%这里为什么不union
  81.         supdofs = union(supdofs1,supdofs2);
  82.     case 'lbracket'
  83.         % L-bracket
  84.         [supnodes,~] = find(X(:,2)==sizey); % Find nodes with y==sizey
  85.         supdofs = union(2*supnodes-1,2*supnodes);
  86.     case 'ubracket'
  87.         % U-bracket
  88.         [supnodes,~] = find(X(:,1)==0); % Find nodes with x==0
  89.         supdofs = union(2*supnodes-1,2*supnodes);
  90. end
  91. alldofs = 1:ndof;
  92. freedofs = setdiff(alldofs,supdofs);
  93. % Null space elimination of supports
  94. N = ones(ndof,1); N(supdofs,1) = 0; Null = spdiags(N,0,ndof,ndof);
  95. %--------------------------------------------------------------
  96. %% PREPARE FILTER
  97. iH = ones(nelem*100,1);
  98. jH = ones(size(iH));
  99. sH = zeros(size(iH));
  100. k = 0;
  101. for e1 = 1:nelem
  102.     x1 = T(e1,7);
  103.     y1 = T(e1,8);
  104.     for e2 = 1:nelem
  105.         x2 = T(e2,7);
  106.         y2 = T(e2,8);
  107.         dist = sqrt((x2-x1)^2+(y2-y1)^2);
  108.         if dist<=rmin
  109.             k = k+1;
  110.             iH(k) = e1;
  111.             jH(k) = e2;
  112.             sH(k) = rmin-dist;
  113.         end
  114.     end
  115. end
  116. H = sparse(iH,jH,sH);
  117. Hs = sum(H,2);
  118. %--------------------------------------------------------------
  119. %% INITIALIZE OPTIMIZATION
  120. maxloop = 120;
  121. volfrac = 0.95;
  122. x = volfrac*ones(nelem,1);
  123. x(T(:,5)==1) = 1e-6; % Voids
  124. x(T(:,5)==2) = 1-1e-6; % Solids
  125. beta = 4;
  126. njumps = 7;
  127. dbeta = (betamaxsimul/beta)^(1/njumps); % Factor for multiplying beta
  128. pace = min(20,ceil(maxloop/(njumps+1)));
  129. loop = 0;
  130. %% INITIALIZE MMA OPTIMIZER
  131. m     = 2;                % The number of general constraints.
  132. n     = nelem;            % The number of design variables x_j.
  133. xmin  = 1e-6*ones(n,1);   % Column vector with the lower bounds for the variables x_j.
  134. xmin(T(:,5)==2) = 1-1e-3; % Lower bound for solids
  135. xmax  = ones(n,1);        % Column vector with the upper bounds for the variables x_j.
  136. xmax(T(:,5)==1) = 1e-3;   % Upper bound for voids
  137. xold1 = x(:);             % xval, one iteration ago (provided that iter>1).
  138. xold2 = x(:);             % xval, two iterations ago (provided that iter>2).
  139. low   = 0*ones(n,1);      % Column vector with the lower asymptotes from the previous iteration (provided that iter>1).
  140. upp   = ones(n,1);        % Column vector with the upper asymptotes from the previous iteration (provided that iter>1).
  141. a0    = 1;                % The constants a_0 in the term a_0*z.
  142. a     = zeros(m,1);       % Column vector with the constants a_i in the terms a_i*z.
  143. c_MMA = 1000*ones(m,1);   % Column vector with the constants c_i in the terms c_i*y_i.
  144. d     = ones(m,1);        % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
  145. %% START ITERATION
  146. Sy = Sytarget;
  147. maxVMscaled = 2*Sytarget;
  148. fval = ones(m,1);
  149. Stats = zeros(maxloop,20);
  150. pcgi1 = 0; pcgi2 = 0;
  151. pcgr1 = 0; pcgr2 = 0;
  152. solver.cgtol11 = solver.cgtol1;
  153. solver.cgtol22 = solver.cgtol2;
  154. %% START ITERATION
  155. while ((loop < maxloop || beta<0.99*betamaxsimul) || (maxVMscaled > Sytarget || fval(1,1) > 0))%fval表示为compliance约束
  156.     loop = loop + 1;
  157.     %% FILTERING OF DESIGN VARIABLES
  158.     if ft == 2
  159.         % Density filter only
  160.         xTilde = (H*x)./Hs;
  161.         xPhys = xTilde;
  162.     elseif ft == 3
  163.         % Density + single projection
  164.         xTilde = (H*x)./Hs;
  165.         xPhys = (tanh(beta*0.5)+tanh(beta*(xTilde-0.5)))/...
  166.             (tanh(beta*0.5)+tanh(beta*(1-0.5)));
  167.     end
  168.     %% FE-ANALYSIS
  169.     sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(Emax-Emin)),64*nelem,1);
  170.     K = sparse(iK,jK,sK);
  171.     K = (K+K')/2;
  172.     % For iterative stress calculation
  173.     stress_data.B = Bmat;
  174.     stress_data.D = Dmat;
  175.     stress_data.X = xPhys;
  176.     stress_data.EDOF = edofMat;
  177.     stress_data.E(1) = Emin;
  178.     stress_data.E(2) = Emax;
  179.     if (solver.type==1)
  180.         % Direct
  181.         U(freedofs) = K(freedofs,freedofs)\F(freedofs);
  182.         flag = 0;
  183.     elseif (solver.type==4)
  184.         % MGCG
  185.         nl = 4;
  186.         Kmg = cell(nl,1);
  187.         Kmg{1,1} = Null'*K*Null - (Null-speye(ndof,ndof));
  188.         for l = 1:nl-1
  189.             Kmg{l+1,1} = Pu{l,1}'*(Kmg{l,1}*Pu{l,1});
  190.         end
  191.         Lfac = chol(Kmg{nl,1},'lower'); Ufac = Lfac';   
  192.         [pcgi1,pcgr1,U] = mgcg_stress(Kmg,F,U,Lfac,Ufac,Pu,nl,1,solver.cgtol11,solver.cgmax1,stress_data);
  193.         flag = 0;
  194.     end
  195.     %%%%%%%%%%%%%%%%
  196.     %% COMPLIANCE AND SENSITIVITY
  197.     ce = sum((U(edofMat)*KE).*U(edofMat),2);
  198.     comp = sum(sum((Emin+xPhys.^penal*(Emax-Emin)).*ce));
  199.     dc = -penal*(Emax-Emin)*xPhys.^(penal-1).*ce;
  200.     dv = ones(nelem,1);
  201.     %% STRESS CONSTRAINT
  202.     % Evaluate stresses
  203.     [appSmax,~,~,~,ADJ,~,dsigvec] = getstress2d(nelem,U,xPhys,edofMat,Emin,Emax,nu,8,helem);
  204.     if (solver.type==1)
  205.         % Direct
  206.         LAM(freedofs) = K(freedofs,freedofs)\ADJ(freedofs);
  207.         pcgi2 = 0; pcgr2 = 0;
  208.     elseif (solver.type==4)
  209.         solver.cgtol22 = solver.cgtol11;
  210.         ADJ = Null*ADJ;
  211.         [pcgi2,pcgr2,LAM] = mgcg_adj(Kmg,ADJ,LAM,Lfac,Ufac,Pu,nl,1,solver.cgtol22,solver.cgmax2,...
  212.             U,edofMat,KE,xPhys,penal,Emax,Emin);
  213.     end
  214.     ce1 = sum(LAM(edofMat)*KE.*U(edofMat),2); % Lam'*K*U
  215.     dsig1 = penal*(Emax-Emin)*xPhys.^(penal-1).*ce1; % Lam'*dKdrho*U
  216.     dsig2 = dsigvec; % dsigaggdrho
  217.     %% FILTERING/MODIFICATION OF SENSITIVITIES
  218.     if  ft == 2
  219.         dv(:) = H*(dv(:)./Hs);
  220.         dsig1(:) = H*(dsig1(:)./Hs);
  221.         dsig2(:) = H*(dsig2(:)./Hs);
  222.         dc(:) = H*(dc(:)./Hs);
  223.     elseif ft == 3
  224.         dxphys = (1 - (tanh(beta*(xTilde(:)-0.5))).^2)*beta / ...
  225.             (tanh(beta*0.5)+tanh(beta*(1-0.5)));
  226.         dv(:) = H*(dv(:).*dxphys(:)./Hs);
  227.         dsig1(:) = H*(dsig1(:).*dxphys(:)./Hs);
  228.         dsig2(:) = H*(dsig2(:).*dxphys(:)./Hs);
  229.         dc(:) = H*(dc(:).*dxphys(:)./Hs);
  230.     end
  231.     dsig = dsig2 + dsig1;
  232.     %% KKT check
  233.     if (loop > 1)
  234.     [~,kktnorm,~] = kktcheck(m,n,xmma,ymma,zmma,lam,...
  235.         xsi,eta,mu_mma,zet,s,xmin,xmax,df0dx,fval,dfdx,a0,a,c_MMA,d);
  236.     else
  237.         kktnorm = 10;
  238.     end
  239.     %% Scaling of constraints (with "dilated" density)
  240.     UxPhys = U;   
  241.     if (1>0)
  242.         % So-called "smooth" scaling
  243.         xDilSharp = (tanh(betamax*0.4)+tanh(betamax*(xTilde-0.4)))/...
  244.             (tanh(betamax*0.4)+tanh(betamax*(1-0.4))); % For measuring "actual" stress
  245.     else
  246.         % So-called "sharp" scaling
  247.         xDilSharp = (tanh(64*0.5)+tanh(64*(xTilde-0.5)))/...
  248.             (tanh(64*0.5)+tanh(64*(1-0.5))); % For measuring "actual" stress
  249.     end
  250.     STRAIN = UxPhys(edofMat)*Bmat';
  251.     STRESS = STRAIN*Dmat;
  252.     Emat = repmat(Emin+xDilSharp(:)*(Emax-Emin),1,3);
  253.     STRESS = STRESS.*Emat;
  254.     VM = STRESS(:,1).^2 - STRESS(:,1).*STRESS(:,2) + STRESS(:,2).^2 + 3*STRESS(:,3).^2;
  255.     VM = VM.^(0.5);
  256.     maxVMscaled = max(VM);
  257.     %% DRAW
  258.     figure(1);
  259.     clf;
  260.     subplot(2,3,1);
  261.     v_img = dsig1;
  262.     top_img = sparse(i_img,j_img,v_img);
  263.     imagesc(top_img);
  264.     axis equal;
  265.     axis tight;
  266.     title('dsig1');
  267.     subplot(2,3,2);
  268.     v_img = dsig2;
  269.     top_img = sparse(i_img,j_img,v_img);
  270.     imagesc(top_img);
  271.     axis equal;
  272.     axis tight;
  273.     title('dsig2');
  274.     subplot(2,3,3);
  275.     v_img = dc;
  276.     top_img = sparse(i_img,j_img,v_img);
  277.     imagesc(top_img);
  278.     axis equal;
  279.     axis tight;
  280.     title('dc');
  281.     subplot(2,3,4);
  282.     v_img = xPhys;
  283.     top_img = sparse(i_img,j_img,v_img);
  284.     imagesc(top_img);
  285.     axis equal;
  286.     axis tight;
  287.     title('xPhys');
  288.     subplot(2,3,5);
  289.     v_img = VM;
  290.     top_img = sparse(i_img,j_img,v_img);
  291.     imagesc(top_img);
  292.     axis equal;
  293.     axis tight;
  294.     title('VM(xDil)');
  295.     drawnow;
  296.     %% METHOD OF MOVING ASYMPTOTES
  297.     xval  = x;
  298.     % minVstCS
  299.     f0val = mean(xPhys);
  300.     if (loop==1)
  301.         scale = 10/f0val;
  302.     end
  303.     df0dx = scale*dv(:)/n;
  304.     if (mod(loop,5) == 1)
  305.         % Scale stress constraint
  306.         ratio = maxVMscaled/Sytarget;
  307.         Sy = appSmax/ratio;
  308.     end
  309.     fval(1,1) = appSmax/Sy - 1;
  310.     dfdx(1,:) = dsig'/Sy;
  311.     if (m > 1) % For testing without compliance
  312.         fval(2,1) = comp/comptarget - 1;
  313.         dfdx(2,:) = dc'/comptarget;
  314.     else
  315.         fval(2,1) = 0;
  316.     end
  317.     fval_mma = fval(1:m,:);
  318.     [xmma,ymma,zmma,lam,xsi,eta,mu_mma,zet,s,low,upp] = ...
  319.         mmasub(m,n,loop,xval,max(xmin,xval-0.1),min(xmax,xval+0.1),xold1,xold2, ...
  320.         f0val,df0dx,fval_mma,dfdx,low,upp,a0,a,c_MMA,d);
  321.     % Update MMA Variables
  322.     xnew     = xmma;
  323.     xold2    = xold1(:);
  324.     xold1    = xval(:);
  325.     change = max(abs(xnew(:)-xval(:)));
  326.     x = xnew;
  327.     %% CONTINUATION
  328.     if (mod(loop,pace) == 0)
  329.         beta = min(beta*dbeta,betamaxsimul);
  330.     end
  331.     %% PRINT RESULTS
  332.     fprintf(' It.:%3i Vol.:%4.3f ch.:%4.3f VMapp:%6.3e VMacc:%6.3e g: %6.3e %6.3e betaHS:%6.3f\n',...
  333.         loop,mean(xPhys(:)),change,appSmax,maxVMscaled,fval',beta);
  334.     fprintf(' State solver %4i %6.3e Adjoint solver %4i %6.3e\n',...
  335.         pcgi1,pcgr1,pcgi2,pcgr2);
  336.     %% KEEP STATS
  337.     Stats(loop,1:15+2*(m-1)) = [f0val fval' appSmax maxVMscaled beta pcgi1 pcgr1 pcgi2 pcgr2...
  338.         ymma' zmma' lam' flag kktnorm];
  339. end
  340. runtime = toc;
  341. save(filename);
  342. end
复制代码


最佳答案


20141303 发表于 6 天前
仅供参考,solver应是结构数组,而type是其一个元素,具体可以查询struct或其他关于结构体的资料
回复此楼

2 条回复


20141303 发表于 6 天前
仅供参考,solver应是结构数组,而type是其一个元素,具体可以查询struct或其他关于结构体的资料
回复此楼

Yintruder 发表于 6 天前
20141303 发表于 2021-5-1 20:12
仅供参考,solver应是结构数组,而type是其一个元素,具体可以查询struct或其他关于结构体的资料 ...

谢谢你!
确实是这样的,他的参数结构数组定义在其他m文件中,我没能理清其中的关系。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
热门教程
站长推荐
快速回复 返回顶部 返回列表