简介这份资源面向拓扑优化初学者与结构优化方向的学生、研究人员提供一套可直接运行的MATLAB入门程序涵盖优化准则法OC与移动渐近线法MMA两种主流求解思路帮助读者理解如何在给定荷载与约束下寻找最优材料分布。压缩包共9个文件全部为m脚本整体约22KB体积轻量却覆盖二维与三维拓扑优化流程包括网格划分、材料属性分配、灵敏度计算与迭代更新等核心环节其中MMA相关子程序可支撑带约束的连续变量优化。目前已有3193人学习下载说明其在入门群体中认可度较高。通过逐行研读并动手修改这些代码读者能掌握OC法与MMA的算法差异、梯度计算与收敛策略并熟悉MATLAB矩阵运算与函数封装为扩展到复杂三维问题或自研优化框架打下基础。1. 从两份不到百行的代码说起拓扑优化入门最该先啃的硬骨头很多人第一次接触结构拓扑优化都是从网上流传的那份“99行程序”开始的。它用不到一百行 MATLAB 代码就把一个连续体结构在给定载荷和边界条件下逐步“掏空”成最优传力路径的全过程跑了出来。配套的还有一份 88 行版本核心区别在于求解器换成了移动渐近线法MMA而 99 行版本用的是优化准则法OC。这两份程序之所以被反复下载不是因为它们工程上多能打而是因为它们把拓扑优化最核心的几件事——设计变量更新、灵敏度分析、密度过滤、收敛判据——全部压缩到了一个能一口气读完的体量里。对于想从理论跨到实现的人来说这比啃完半本专著再动手要高效得多。它适合已经学过有限元基础、想亲手改参数看结构变化的人也适合需要快速验证一个新想法、不想从零搭框架的研究者。2. OC 与 MMA 的算法骨架为什么 99 行和 88 行差在这2.1 优化准则法OC的更新逻辑与适用边界OC 的核心思想很直接在体积约束下让每个单元的材料分布满足最优性条件。具体到代码里它构造一个拉格朗日函数对设计变量求导后令其满足 KKT 条件最终推导出一个乘子式的更新公式。99 行程序里那段xnew max(0, max(x-move, min(1, min(xmove, x.*sqrt(-dc./lambda)))))就是 OC 的典型实现。其中dc是目标函数对设计变量的灵敏度lambda是通过二分法找出的拉格朗日乘子move是每步允许的最大变化量通常取 0.2。这个更新式的物理含义是灵敏度越大的单元即增加它对刚度贡献越大越倾向于保留材料灵敏度小甚至为负的单元则被削减。OC 的优点是收敛快、每步计算量小尤其适合单约束体积约束问题。但它的边界也很明显只能处理单约束且更新公式的推导依赖于特定形式的优化问题。一旦你加入多个约束比如同时限制体积和某个局部应力OC 的乘子法框架就需要大改甚至直接失效。很多新手改代码时想加个位移约束结果发现 OC 那套二分法完全对不上号这就是没搞清楚适用边界。2.2 MMA 的渐近线机制与多约束扩展性88 行版本换成了 MMA代码量反而少了因为 MMA 把更多逻辑封装进了子问题求解。MMA 的核心是在每步迭代中构造一个近似子问题用渐近线参数L和U控制近似函数的曲率然后求解这个显式的凸子问题得到更新方向。代码里你会看到low和upp两个数组它们就是渐近线每次迭代根据设计变量的变化幅度自适应调整。如果某步更新后变量撞到了渐近线边界下次迭代就把渐近线往外扩如果变化平缓就收紧。MMA 最大的价值在于它能统一处理多约束问题。你只需要把约束的灵敏度和约束值传进去子问题求解器会自动处理。这也是为什么工业级拓扑优化软件几乎都用 MMA 或它的变体如 GCMMA。但代价是每步计算量比 OC 大且渐近线参数需要调——asyinit、asyincr、asydecr这几个参数设不好收敛曲线会剧烈震荡甚至发散。我一般建议新手先用默认值跑通再逐步调asyinit从 0.5 降到 0.3 观察收敛速度变化。2.3 两种算法的代码结构对比把两份代码并排看结构差异一目了然。99 行版本的循环体里有限元分析、灵敏度计算、过滤、OC 更新、收敛判断全部平铺88 行版本则把更新步骤替换为对mmasub的调用主循环更干净。下面这段是 OC 更新的核心逻辑我加了注释方便对照% OC 更新核心段99行版本节选 l1 0; l2 1e9; % 拉格朗日乘子的上下界 while (l2-l1)/(l1l2) 1e-3 % 二分法找乘子 lmid 0.5*(l2l1); xnew max(0, max(x-move, min(1, min(xmove, x.*sqrt(-dc./lmid))))); if sum(sum(xnew)) - volfrac*nelx*nely 0 l1 lmid; % 体积超了乘子取大 else l2 lmid; % 体积没超乘子取小 end end这段代码里volfrac是目标体积分数nelx和nely是单元划分数量。二分法的收敛容差1e-3控制乘子精度太松会导致体积约束漂移太紧则增加内层循环次数。实际跑的时候如果发现最终体积分数总是比设定值高一点先把容差调到1e-4试试。MMA 那边对应的调用则是% MMA 子问题调用88行版本节选 [xnew, ynew, znew, lamnew] mmasub(m, n, iter, xval, xmin, xmax, ... xold1, xold2, f0val, df0dx, fval, dfdx, low, upp, a0, a, c, d);参数m是约束个数n是设计变量个数low和upp是渐近线。注意xmin和xmax通常设为 0.001 和 1而不是 0 和 1——因为灵敏度计算里会出现除以密度的情况密度为 0 会导致数值问题。这个 0.001 的下界就是所谓的“空材料”近似是 SIMP 模型的标准做法。3. 跑通第一张拓扑图从参数设置到结果解读3.1 有限元网格与边界条件的代码映射拿到代码后第一件事是改边界条件。以经典的悬臂梁为例99 行程序里载荷施加在右下角节点左边界全约束。代码里对应的是% 边界条件与载荷以悬臂梁为例 ndof 2*(nelx1)*(nely1); % 总自由度 fixeddofs 1:2*(nely1); % 左边界全约束 alldofs 1:ndof; freedofs setdiff(alldofs, fixeddofs); % 载荷施加在右下角 loaddof 2*(nelx1)*(nely1) - 1; F sparse(loaddof, 1, 1, ndof, 1);这里fixeddofs的编号规则是节点按列优先排列每个节点两个自由度x 和 y 方向。左边界节点编号从 1 到nely1对应的 x 方向自由度是1, 3, 5, ...y 方向是2, 4, 6, ...。所以1:2*(nely1)正好覆盖左边界所有自由度。如果你要改成两端固支就把右边界也加进去但注意载荷位置要相应调整否则结构会变成纯受压状态拓扑结果会是一堆无意义的材料堆积。3.2 密度过滤与灵敏度修正直接对单元密度做优化会出现棋盘格现象——相邻单元密度在 0 和 1 之间高频振荡。99 行和 88 行程序都用了灵敏度过滤来解决核心代码是% 灵敏度过滤99行版本节选 ft 2; % 过滤类型1为灵敏度过滤2为密度过滤 rmin 2.5; % 过滤半径 if ft 1 dc conv2(dc.*x, ones(nely,nelx)/(nely*nelx), same) ... ./ (x.*conv2(ones(nely,nelx), ones(nely,nelx)/(nely*nelx), same)); endrmin是最关键的参数。设得太小棋盘格压不住设得太大结构会过度平滑细杆件消失。对于 60x30 的网格rmin2.5是个稳妥的起点。如果你把网格加密到 120x60rmin也要相应增大到 4 左右否则过滤效果会被稀释。注意conv2的边界处理用的是same这会导致边界附近的过滤权重偏小实际有效过滤半径比设定值略小这是常见做法不用刻意修正。3.3 收敛判据与迭代终止程序默认的收敛判据是设计变量最大变化量小于 0.01change max(abs(xnew(:)-x(:))); if change 0.01 break; end这个阈值对大多数二维问题够用但如果你发现结构在 80 步左右还在缓慢变化别急着停——可能是move设得太小导致每步走得太保守。把move从 0.2 提到 0.3 试试收敛会快一些但太大又容易震荡。我一般会同时监控目标函数柔度的变化如果连续 5 步柔度下降不到 0.1%即使change没到阈值也可以认为收敛了。另外注意OC 的二分法内循环如果没收敛外层迭代也会受影响表现为体积分数在目标值附近来回跳这时候先检查二分容差和move的配合。4. 避坑与排查那些让程序跑不出结果的细节4.1 现象拓扑结果全是灰色没有黑白分明原因通常是过滤半径过大或迭代步数不够。rmin设到 5 以上密度场会被过度平滑所有单元都落在 0.5 附近。解决方法是把rmin降到 2 左右同时检查ft是否设成了 1灵敏度过滤而不是 2密度过滤。密度过滤在 99 行版本里需要额外处理如果直接套用灵敏度过滤的代码结果会偏灰。4.2 现象体积分数始终不收敛在目标值上下大幅波动先看 OC 的二分法容差是不是太松。把1e-3改成1e-4能改善但更常见的原因是move太大。move0.5时每步允许密度变化一半体积约束很容易被冲过头。改回 0.2 再跑。如果是 MMA检查asyinit是不是设得太小渐近线收得太紧会导致子问题近似失真更新方向乱跳。4.3 现象柔度曲线先降后升结构越跑越差这是典型的灵敏度符号错误或过滤方向反了。检查dc的符号对于柔度最小化问题dc应该是负的因为增加材料会降低柔度OC 更新式里用的是-dc。如果你在改代码时不小心把符号弄反结构会朝着增加柔度的方向演化。另一个可能是过滤时把dc.*x写成了dc./x这会让低密度区域的灵敏度被异常放大。4.4 现象程序报错“矩阵维度不匹配”九成是网格尺寸和载荷位置对不上。nelx和nely改了之后loaddof的计算必须同步更新。比如你把nelx从 60 改成 80右下角节点的编号就变了原来的loaddof公式如果写死了具体数值就会出错。建议始终用2*(nelx1)*(nely1)-1这种表达式而不是硬编码。另外fixeddofs的范围也要跟着nely走别漏改。4.5 现象MMA 版本跑几步就提示“子问题无解”先检查xmin是不是设成了 0。MMA 的子问题求解器要求设计变量有严格正的下界因为渐近线更新公式里会出现除以x-xmin的项。把xmin改成1e-3或0.001。如果还不行看a0和a这些约束参数有没有传错——a0是目标函数的线性化系数通常设为 0 或 1具体取决于你的问题形式。88 行版本里这些参数在循环外初始化改约束时别忘了同步改。5. 从跑通到改出自己需要的东西几个能直接抄的进阶操作跑通默认算例只是起点真正要用起来得知道怎么把问题改成自己的。第一个操作是换目标函数。默认是最小化柔度如果你想做 compliant mechanism 设计把目标改成输出位移最大化只需要把灵敏度dc的符号反过来同时把c的计算从sum(sum(E0.*x.^penal.*ue.*ue))改成对应的输出位移表达式。注意这时候体积约束依然保留但 OC 的更新公式需要重新推导——这也是为什么做机构设计时建议直接用 MMA 版本省去推导麻烦。第二个操作是加被动单元。比如你想固定某些区域必须保留材料在每次更新后强制x(passive) 1同时把这些单元的灵敏度置零防止它们被优化器改动。代码里加两行% 被动单元处理在更新后、过滤前插入 passive false(nely, nelx); passive(1:2, :) true; % 顶部两行必须保留 x(passive) 1; dc(passive) 0;第三个操作是导出结果。MATLAB 里用imagesc看结果方便但要导入其他软件做后续分析得输出 STL 或网格文件。一个偷懒的办法是把密度矩阵存成 CSV再用 Python 的numpy和matplotlib重建等值面。我常用的是writematrix(x, density.csv)然后在 Python 里用plt.contourf画出来阈值取 0.5。如果要更精确的边界用skimage.measure.marching_squares提取轮廓再导出 DXF。最后一个技巧是参数扫描。与其手动改一个跑一次不如写个外层循环批量跑不同volfrac和rmin组合把最终柔度和迭代步数存下来。这样你能快速看出哪个参数区间结构最合理。我一般会固定rmin2.5让volfrac从 0.3 到 0.6 步进 0.1 跑五组对比柔度下降幅度和结构清晰度。从那以后我每次改代码加新约束都强制先用小网格30x15跑一遍确认逻辑没问题再放大网格——直接上大网格一旦出错等半天才看到结果后悔药都没得吃。希望帮到你。本文还有配套的精品资源点击获取
简介这份资源面向拓扑优化初学者与结构优化方向的学生、研究人员提供一套可直接运行的MATLAB入门程序涵盖优化准则法OC与移动渐近线法MMA两种主流求解思路帮助读者理解如何在给定荷载与约束下寻找最优材料分布。压缩包共9个文件全部为m脚本整体约22KB体积轻量却覆盖二维与三维拓扑优化流程包括网格划分、材料属性分配、灵敏度计算与迭代更新等核心环节其中MMA相关子程序可支撑带约束的连续变量优化。目前已有3193人学习下载说明其在入门群体中认可度较高。通过逐行研读并动手修改这些代码读者能掌握OC法与MMA的算法差异、梯度计算与收敛策略并熟悉MATLAB矩阵运算与函数封装为扩展到复杂三维问题或自研优化框架打下基础。1. 从两份不到百行的代码说起拓扑优化入门最该先啃的硬骨头很多人第一次接触结构拓扑优化都是从网上流传的那份“99行程序”开始的。它用不到一百行 MATLAB 代码就把一个连续体结构在给定载荷和边界条件下逐步“掏空”成最优传力路径的全过程跑了出来。配套的还有一份 88 行版本核心区别在于求解器换成了移动渐近线法MMA而 99 行版本用的是优化准则法OC。这两份程序之所以被反复下载不是因为它们工程上多能打而是因为它们把拓扑优化最核心的几件事——设计变量更新、灵敏度分析、密度过滤、收敛判据——全部压缩到了一个能一口气读完的体量里。对于想从理论跨到实现的人来说这比啃完半本专著再动手要高效得多。它适合已经学过有限元基础、想亲手改参数看结构变化的人也适合需要快速验证一个新想法、不想从零搭框架的研究者。2. OC 与 MMA 的算法骨架为什么 99 行和 88 行差在这2.1 优化准则法OC的更新逻辑与适用边界OC 的核心思想很直接在体积约束下让每个单元的材料分布满足最优性条件。具体到代码里它构造一个拉格朗日函数对设计变量求导后令其满足 KKT 条件最终推导出一个乘子式的更新公式。99 行程序里那段xnew max(0, max(x-move, min(1, min(xmove, x.*sqrt(-dc./lambda)))))就是 OC 的典型实现。其中dc是目标函数对设计变量的灵敏度lambda是通过二分法找出的拉格朗日乘子move是每步允许的最大变化量通常取 0.2。这个更新式的物理含义是灵敏度越大的单元即增加它对刚度贡献越大越倾向于保留材料灵敏度小甚至为负的单元则被削减。OC 的优点是收敛快、每步计算量小尤其适合单约束体积约束问题。但它的边界也很明显只能处理单约束且更新公式的推导依赖于特定形式的优化问题。一旦你加入多个约束比如同时限制体积和某个局部应力OC 的乘子法框架就需要大改甚至直接失效。很多新手改代码时想加个位移约束结果发现 OC 那套二分法完全对不上号这就是没搞清楚适用边界。2.2 MMA 的渐近线机制与多约束扩展性88 行版本换成了 MMA代码量反而少了因为 MMA 把更多逻辑封装进了子问题求解。MMA 的核心是在每步迭代中构造一个近似子问题用渐近线参数L和U控制近似函数的曲率然后求解这个显式的凸子问题得到更新方向。代码里你会看到low和upp两个数组它们就是渐近线每次迭代根据设计变量的变化幅度自适应调整。如果某步更新后变量撞到了渐近线边界下次迭代就把渐近线往外扩如果变化平缓就收紧。MMA 最大的价值在于它能统一处理多约束问题。你只需要把约束的灵敏度和约束值传进去子问题求解器会自动处理。这也是为什么工业级拓扑优化软件几乎都用 MMA 或它的变体如 GCMMA。但代价是每步计算量比 OC 大且渐近线参数需要调——asyinit、asyincr、asydecr这几个参数设不好收敛曲线会剧烈震荡甚至发散。我一般建议新手先用默认值跑通再逐步调asyinit从 0.5 降到 0.3 观察收敛速度变化。2.3 两种算法的代码结构对比把两份代码并排看结构差异一目了然。99 行版本的循环体里有限元分析、灵敏度计算、过滤、OC 更新、收敛判断全部平铺88 行版本则把更新步骤替换为对mmasub的调用主循环更干净。下面这段是 OC 更新的核心逻辑我加了注释方便对照% OC 更新核心段99行版本节选 l1 0; l2 1e9; % 拉格朗日乘子的上下界 while (l2-l1)/(l1l2) 1e-3 % 二分法找乘子 lmid 0.5*(l2l1); xnew max(0, max(x-move, min(1, min(xmove, x.*sqrt(-dc./lmid))))); if sum(sum(xnew)) - volfrac*nelx*nely 0 l1 lmid; % 体积超了乘子取大 else l2 lmid; % 体积没超乘子取小 end end这段代码里volfrac是目标体积分数nelx和nely是单元划分数量。二分法的收敛容差1e-3控制乘子精度太松会导致体积约束漂移太紧则增加内层循环次数。实际跑的时候如果发现最终体积分数总是比设定值高一点先把容差调到1e-4试试。MMA 那边对应的调用则是% MMA 子问题调用88行版本节选 [xnew, ynew, znew, lamnew] mmasub(m, n, iter, xval, xmin, xmax, ... xold1, xold2, f0val, df0dx, fval, dfdx, low, upp, a0, a, c, d);参数m是约束个数n是设计变量个数low和upp是渐近线。注意xmin和xmax通常设为 0.001 和 1而不是 0 和 1——因为灵敏度计算里会出现除以密度的情况密度为 0 会导致数值问题。这个 0.001 的下界就是所谓的“空材料”近似是 SIMP 模型的标准做法。3. 跑通第一张拓扑图从参数设置到结果解读3.1 有限元网格与边界条件的代码映射拿到代码后第一件事是改边界条件。以经典的悬臂梁为例99 行程序里载荷施加在右下角节点左边界全约束。代码里对应的是% 边界条件与载荷以悬臂梁为例 ndof 2*(nelx1)*(nely1); % 总自由度 fixeddofs 1:2*(nely1); % 左边界全约束 alldofs 1:ndof; freedofs setdiff(alldofs, fixeddofs); % 载荷施加在右下角 loaddof 2*(nelx1)*(nely1) - 1; F sparse(loaddof, 1, 1, ndof, 1);这里fixeddofs的编号规则是节点按列优先排列每个节点两个自由度x 和 y 方向。左边界节点编号从 1 到nely1对应的 x 方向自由度是1, 3, 5, ...y 方向是2, 4, 6, ...。所以1:2*(nely1)正好覆盖左边界所有自由度。如果你要改成两端固支就把右边界也加进去但注意载荷位置要相应调整否则结构会变成纯受压状态拓扑结果会是一堆无意义的材料堆积。3.2 密度过滤与灵敏度修正直接对单元密度做优化会出现棋盘格现象——相邻单元密度在 0 和 1 之间高频振荡。99 行和 88 行程序都用了灵敏度过滤来解决核心代码是% 灵敏度过滤99行版本节选 ft 2; % 过滤类型1为灵敏度过滤2为密度过滤 rmin 2.5; % 过滤半径 if ft 1 dc conv2(dc.*x, ones(nely,nelx)/(nely*nelx), same) ... ./ (x.*conv2(ones(nely,nelx), ones(nely,nelx)/(nely*nelx), same)); endrmin是最关键的参数。设得太小棋盘格压不住设得太大结构会过度平滑细杆件消失。对于 60x30 的网格rmin2.5是个稳妥的起点。如果你把网格加密到 120x60rmin也要相应增大到 4 左右否则过滤效果会被稀释。注意conv2的边界处理用的是same这会导致边界附近的过滤权重偏小实际有效过滤半径比设定值略小这是常见做法不用刻意修正。3.3 收敛判据与迭代终止程序默认的收敛判据是设计变量最大变化量小于 0.01change max(abs(xnew(:)-x(:))); if change 0.01 break; end这个阈值对大多数二维问题够用但如果你发现结构在 80 步左右还在缓慢变化别急着停——可能是move设得太小导致每步走得太保守。把move从 0.2 提到 0.3 试试收敛会快一些但太大又容易震荡。我一般会同时监控目标函数柔度的变化如果连续 5 步柔度下降不到 0.1%即使change没到阈值也可以认为收敛了。另外注意OC 的二分法内循环如果没收敛外层迭代也会受影响表现为体积分数在目标值附近来回跳这时候先检查二分容差和move的配合。4. 避坑与排查那些让程序跑不出结果的细节4.1 现象拓扑结果全是灰色没有黑白分明原因通常是过滤半径过大或迭代步数不够。rmin设到 5 以上密度场会被过度平滑所有单元都落在 0.5 附近。解决方法是把rmin降到 2 左右同时检查ft是否设成了 1灵敏度过滤而不是 2密度过滤。密度过滤在 99 行版本里需要额外处理如果直接套用灵敏度过滤的代码结果会偏灰。4.2 现象体积分数始终不收敛在目标值上下大幅波动先看 OC 的二分法容差是不是太松。把1e-3改成1e-4能改善但更常见的原因是move太大。move0.5时每步允许密度变化一半体积约束很容易被冲过头。改回 0.2 再跑。如果是 MMA检查asyinit是不是设得太小渐近线收得太紧会导致子问题近似失真更新方向乱跳。4.3 现象柔度曲线先降后升结构越跑越差这是典型的灵敏度符号错误或过滤方向反了。检查dc的符号对于柔度最小化问题dc应该是负的因为增加材料会降低柔度OC 更新式里用的是-dc。如果你在改代码时不小心把符号弄反结构会朝着增加柔度的方向演化。另一个可能是过滤时把dc.*x写成了dc./x这会让低密度区域的灵敏度被异常放大。4.4 现象程序报错“矩阵维度不匹配”九成是网格尺寸和载荷位置对不上。nelx和nely改了之后loaddof的计算必须同步更新。比如你把nelx从 60 改成 80右下角节点的编号就变了原来的loaddof公式如果写死了具体数值就会出错。建议始终用2*(nelx1)*(nely1)-1这种表达式而不是硬编码。另外fixeddofs的范围也要跟着nely走别漏改。4.5 现象MMA 版本跑几步就提示“子问题无解”先检查xmin是不是设成了 0。MMA 的子问题求解器要求设计变量有严格正的下界因为渐近线更新公式里会出现除以x-xmin的项。把xmin改成1e-3或0.001。如果还不行看a0和a这些约束参数有没有传错——a0是目标函数的线性化系数通常设为 0 或 1具体取决于你的问题形式。88 行版本里这些参数在循环外初始化改约束时别忘了同步改。5. 从跑通到改出自己需要的东西几个能直接抄的进阶操作跑通默认算例只是起点真正要用起来得知道怎么把问题改成自己的。第一个操作是换目标函数。默认是最小化柔度如果你想做 compliant mechanism 设计把目标改成输出位移最大化只需要把灵敏度dc的符号反过来同时把c的计算从sum(sum(E0.*x.^penal.*ue.*ue))改成对应的输出位移表达式。注意这时候体积约束依然保留但 OC 的更新公式需要重新推导——这也是为什么做机构设计时建议直接用 MMA 版本省去推导麻烦。第二个操作是加被动单元。比如你想固定某些区域必须保留材料在每次更新后强制x(passive) 1同时把这些单元的灵敏度置零防止它们被优化器改动。代码里加两行% 被动单元处理在更新后、过滤前插入 passive false(nely, nelx); passive(1:2, :) true; % 顶部两行必须保留 x(passive) 1; dc(passive) 0;第三个操作是导出结果。MATLAB 里用imagesc看结果方便但要导入其他软件做后续分析得输出 STL 或网格文件。一个偷懒的办法是把密度矩阵存成 CSV再用 Python 的numpy和matplotlib重建等值面。我常用的是writematrix(x, density.csv)然后在 Python 里用plt.contourf画出来阈值取 0.5。如果要更精确的边界用skimage.measure.marching_squares提取轮廓再导出 DXF。最后一个技巧是参数扫描。与其手动改一个跑一次不如写个外层循环批量跑不同volfrac和rmin组合把最终柔度和迭代步数存下来。这样你能快速看出哪个参数区间结构最合理。我一般会固定rmin2.5让volfrac从 0.3 到 0.6 步进 0.1 跑五组对比柔度下降幅度和结构清晰度。从那以后我每次改代码加新约束都强制先用小网格30x15跑一遍确认逻辑没问题再放大网格——直接上大网格一旦出错等半天才看到结果后悔药都没得吃。希望帮到你。本文还有配套的精品资源点击获取
