IGA拓扑优化
一、概述
传统拓扑优化（如基于有限元方法）在结构建模与优化间存在割裂问题，几何建模依赖 CAD，而分析依赖 FEM，这种 几何与分析的不一致性 会引发误差积累与建模复杂性。为解决这一问题，论文引入 等几何分析（IGA），将 CAD 中常用的 NURBS（非均匀有理 B 样条）曲线直接用于分析建模，实现建模与分析的统一。
本论文基于 MATLAB 编写了一个轻量、紧凑、高效的 2D 等几何拓扑优化工具 IgaTop
二、模块介绍
| 模块名 | 功能说明 |
|---|---|
Pre_IGA |
控制点序列构造、自由度编号 |
Geom_Mod |
结构几何建模，生成 NURBS 模型 |
Boun_Cond |
设置 Dirichlet 和 Neumann 边界 |
Assemble |
组装全局刚度矩阵 |
FEA |
结构有限元求解 |
Sensitivity |
灵敏度分析 |
Filter |
应用灵敏度过滤 |
Optimization |
使用 MMA（移动渐近线）优化迭代 |
三、论文解读
\[
\text{Find: } \rho_{i,j} \quad (i=1,2,\dots,n;\ j=1,2,\dots,m)
\]
\[
\text{Min: } J(\mathbf{u}, \mathcal{X}) = \frac{1}{2} \int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{u})^T\, \mathbf{D}(\mathcal{X}(\xi, \eta))\, \boldsymbol{\varepsilon}(\mathbf{u})\, d\Omega
\]
\[
\text{s.t. }
\begin{cases}
G(\mathcal{X}) = \frac{1}{|\Omega|} \int_\Omega \mathcal{X}(\xi, \eta)\, v_0\, d\Omega - V_{\text{max}} \leq 0 & \text{（体积约束）} \\
a(\mathbf{u}, \delta \mathbf{u}) = l(\delta \mathbf{u}),\quad \mathbf{u}|_{\Gamma_D} = \mathbf{g},\ \forall\, \delta \mathbf{u} \in H^1(\Omega) & \text{（平衡方程）} \\
0 \leq \rho_{i,j} \leq 1 & \text{（设计变量约束）}
\end{cases}
\]
以上为论文中提到的ITO的数学表示模型，表达一个典型的基于密度的最小柔度问题 其中
- \(ρ_{i,j}\)表示设计变量
- \(J(\mathbf{u}, \mathcal{X}) = \dfrac{1}{2} \int_\Omega \boldsymbol{\varepsilon}(\mathbf{u})^T \mathbf{D}(\mathcal{X}(\xi,\eta)) \boldsymbol{\varepsilon}(\mathbf{u})\, d\Omega \\\\[1.2ex]\)表示目标函数
- u：位移场 \(\epsilon\)：应变
- D：弹性刚度矩阵，与密度函数 \(\mathcal{X}\)有关
- \(G(\mathcal{X}) = \frac{1}{|\Omega|} \int_\Omega \mathcal{X}(\xi,\eta)\, v_0 \, d\Omega - V_{\text{max}} \leq 0\)
- 表示体积约束，限制材料总体积不超过某个阈值\(V_{max}\)
- \(a(\mathbf{u}, \delta \mathbf{u}) = l(\delta \mathbf{u}),\ \mathbf{u}|_{\Gamma_D} = \mathbf{g},\ \forall\, \delta \mathbf{u} \in H^1(\Omega)\)
- 对应结构的弱形式平衡方程（虚功原理）
- \(\Gamma_D\)：Dirichlet 边界，施加位移条件
- 最后一个为设计变量约束
母空间
母空间是指一个标准区间，例如：
- 1D：[-1,1]
- 2D：[-1,1] \(\times\) [-1,1]
参数空间
定义Nurbs基函数变量的标准空间
物理空间
将 NURBS 基函数定义的“参数空间”映射到结构所在的实际物理空间。使用 NURBS 基函数对参数空间点进行加权组合，实现曲线/面的光滑表达：\(\mathbf{x}(\xi, \eta) = \sum_i R_i(\xi, \eta)\, \mathbf{P}_i\)
在计算刚度矩阵时，需要使用到高斯积分，这被定义到母空间
IGA需要使用两个映射：
母空间 → 参数空间 → 物理空间
上图为使用 IGA 实现拓扑优化的整体框架流程图
3.1 Geom_Mod
在IGA中，二维Nurbs曲面使用如下公式表示
\[
\mathbf{x}(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} R_{i,j}(\xi, \eta)\, \mathbf{P}_{i,j}
\]
其中
- \(R_{i,j}(\xi,\eta)\)：Nurbs基函数
- \(P_{i,j}\)：控制点坐标
Nurbs基函数如下：
\[
R_{i,j}(\xi, \eta) = \frac{N_{i,p}(\xi)\, M_{j,q}(\eta)\, w_{i,j}}{\sum_{k=1}^{n} \sum_{l=1}^{m} N_{k,p}(\xi)\, M_{l,q}(\eta)\, w_{k,l}}
\]
- \(N_{i,p}、M_{j,q}\)：两个方向的B样条基函数
- \(w_{i,j}\)：控制点权重
对于0阶B样条基函数，有如下性质
\[
N_{i,0}(\xi) =
\begin{cases}
1, & \text{if } \xi_i \leq \xi < \xi_{i+1} \\
0, & \text{otherwise}
\end{cases}
\]
对于高阶B样条基函数，遵从下面的递归公式
\[
N_{i,p}(\xi) = \frac{\xi - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\xi)
+ \frac{\xi_{i+p+1} - \xi}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\xi)
\]
NURBS工具箱:
- nrbmak：基于初始结点向量和控制点构建NURBS曲面
- nrbdegelev：提升NURBS基函数在第二个参数方向上的阶数
- nrbkntins： 在初始结点向量中均匀插入一系列新结点
- nrbbasisfun：用于在指定参数位置计算非零基函数和其控制点索引。
- nrbeval: 通过参数坐标计算实际物理坐标
3.2 Pre_IGA
根据Nurbs曲面得出返回三个参数——[CtrPts,Ele,GauPts]
CtrPts：控制点相关信息
| 字段名 | 中文解释 |
|---|---|
CtrPts.Cordis |
控制点在物理空间中的笛卡尔坐标，表示为 \((x, y, z, \omega)\) |
CtrPts.Num |
控制点的总数 |
CtrPts.NumU |
第一参数方向上的控制点数量 |
CtrPts.NumV |
第二参数方向上的控制点数量 |
CtrPts.Seque |
所有控制点的编号矩阵 |
Ele：IGA的单元
| 编号 | 字段名 | 含义 |
|---|---|---|
| 1 | NumU |
第一参数方向（如 ξ）上 IGA 元素的总数（元素片段个数） |
| 2 | NumV |
第二参数方向（如 η）上 IGA 元素的总数 |
| 3 | Num |
IGA 总元素数（= NumU × NumV） |
| 4 | Seque |
所有 IGA 元素的编号（按“左到右、下到上”的顺序） |
| 5 | KnotsU |
每个 IGA 元素在 第一参数方向上的 knot 跨度（区间） |
| 6 | KnotsV |
每个 IGA 元素在 第二参数方向上的 knot 跨度（区间） |
| 7 | CtrPtsNum |
每个元素中对其有影响的控制点数 = 非零 NURBS 基函数阶数（次数+1） |
| 8 | CtrPtsNumU |
每个元素在第一参数方向上的控制点数 |
| 9 | CtrPtsNumV |
每个元素在第二参数方向上的控制点数 |
| 10 | CtrPtsCon |
控制点连接矩阵：第 i 行列出第 i 个元素的控制点编号（共 CtrPtsNum 个）类似 FEM 中的 edofMat |
| 11 | GauPtsNum |
每个元素的 高斯积分点数量，用于数值积分 |
Warning
其中Seque的顺序和控制点编号的顺序一致
GauPts：高斯积分点相关信息
| 字段名 | 中文解释 |
|---|---|
QuaPts |
母空间中高斯积分点的坐标\((\hat{\xi}, \hat{\eta}\)） |
Weigh |
每个高斯点的积分权重 |
Num |
所有高斯积分点的总数 |
CorU |
高斯点从母空间映射到参数空间时第一方向的 knot 对应 |
CorV |
高斯点从母空间映射到参数空间时第二方向的 knot 对应 |
Seque |
每个 IGA 元素对应的高斯点编号矩阵 |
Warning
高斯积分点被定义在母空间中，需要转换到参数空间，转换如下
\(CorU = ((u2 - u1) * \hat{\xi} + (u2 + u1)) / 2\)
\(CorV = ((v2 - v1) * \hat{\eta} + (v2 + v1)) / 2\)
高斯积分点的信息通过Guadrature函数给出
1 | |
其中
- quadorder：高斯积分的阶数（例如 2阶、3阶）
- dim: 积分的维度，当前只处理了 dim = 2
- quadpoint: 高斯积分点坐标（二维）
- quadweight: 每个高斯点的积分权重
3.3 Boun_Cond
用于定义狄利克雷边界和纽曼边界条件
- 狄利克雷边界：固定
- 纽曼边界：载荷
1 | |
| 参数名 | 含义 |
|---|---|
CtrPts |
控制点信息结构体，来自 Pre_IGA，包含控制点坐标 Cordis、数量 Num、顺序编号 Seque 等，用于定义几何边界和施加载荷位置。 |
BoundCon |
整型编号，用于表示使用哪一个边界条件案例（共支持 5 种，如悬臂梁、MBB 梁、L 梁等），影响 Dirichlet/Neumann 条件的设置。 |
NURBS |
NURBS 曲面结构体，包含 knot 矢量、阶数、控制点、权重等信息，是计算基函数值与导数的关键。 |
Dofs.Num |
总自由度数量，等于控制点数 × 每个控制点的自由度（如二维问题中为 2），用于初始化力向量 F 的维度。 |
3.4 DDF
本章用于初始化控制密度以及获取高斯积分点处的DDF(密度分布函数)
- 分配一系列定义在控制点处的离散密度，称为控制密度
- 定义平滑机制以提高控制密度的平滑度
- 将NURBS基函数（用于构建初始结构几何的NURBS曲面）与平滑后的控制密度进行线性组合，从而为整个结构几何开发相应的密度分布函数(DDF)。
其数学表达式为：
\[
\mathcal{X}(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} R_{i,j}^{p,q}(\xi, \eta)\, \tilde{\rho}_{i,j}
\]
DDF：简单点说就是通过Nurbs基函数对控制点密度进行加权累加，得到一个光滑的密度系统，每个点的密度都可以用对应控制点密度获得
1 2 3 4 5 6 7 8 9 10 11 | |
3.5 Shep_Fun
定义平滑机制
1 | |
- 输入控制点和影响半径
- 返回加权矩阵Sh和每个控制点的加权和Hs
定义Shepard平滑函数如下：
\[
\tilde{\rho}_{i,j} =
\sum_{i=1}^{\mathcal{N}} \sum_{j=1}^{\mathcal{M}} \psi(\rho_{i,j}) \rho_{i,j} =
\sum_{i=1}^{\mathcal{N}} \sum_{j=1}^{\mathcal{M}}
\left( \frac{w(\rho_{i,j})}
{\sum_{\hat{i}=1}^{\mathcal{N}} \sum_{\hat{j}=1}^{\mathcal{M}} w(\rho_{\hat{i}, \hat{j}})} \right)
\rho_{i,j}=\frac{ \sum_{i=1}^{N} \sum_{j=1}^{M} w(\rho_{i,j}) \cdot \rho_{i,j} }
{ \sum_{i=1}^{N} \sum_{j=1}^{M} w(\rho_{i,j}) }
\]
\(\psi(\rho_{i,j})=\frac{w(\rho_{i,j})} {\sum_{\hat{i}=1}^{\mathcal{N}} \sum_{\hat{j}=1}^{\mathcal{M}} w(\rho_{\hat{i}, \hat{j}})}\) 即为平滑函数
- \(\rho_{i,j}\)：原始控制密度
- \(\tilde{\rho}_{i,j}\)：平滑后的控制密度；
- \(w(\rho_{i,j})\)：权重函数，反映距离和影响范围
其中
\[w(v)=(1-v)^6+（35v^2+18v+3）\\ v = \frac{r}{r_{m}}
\]
- \(r_m\)：为影响区域的半径。
- \(r\)：当前控制密度与局部支撑域中其他控制密度之间的欧几里得距离
如下图所示
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | |
3.6 structural responses
本章用于求解结构响应，分为三个部分：Stiff_Ele2D、Stiff_Ass2D、Solving
- Stiff_Ele2D：计算所有等几何分析单元的刚度矩阵
- Stiff_Ass2D：组装所有等几何分析单元刚度矩阵
- Solving：解决结构响应问题
Stiff_Ele2D：
等几何分析单元刚度矩阵计算如下
\[
K_e = \int_{\Omega_e} B^T D B \, d\Omega_e
\]
需要使用高斯积分计算，由于高斯积分被定义在母空间中，故需要进行母空间->参数空间的映射，通过构建雅可比矩阵实现，如下式：
\[
K_e = \int_{\tilde{\Omega}_e} B^T D B |J_1||J_2| \, d\tilde{\Omega}_e
\]
\(|J_1|、|J_2|\)分别是x和y方向的雅可比矩阵
1 | |
| 参数 | 含义 |
|---|---|
X |
结构体数组，包含设计变量密度的信息。具体包含： - CtrPts（控制点密度）- GauPts（高斯积分点密度）- DDF（设计域中的密度场） |
penal |
惩罚参数，SIMP模型中通常取3，用于控制材料插值的非线性。 |
Emin |
空材料的最小杨氏模量，避免数值奇异性。 |
DH |
材料本构矩阵D（如平面应力或应变条件下的弹性矩阵）。 |
CtrPts |
控制点信息（NURBS控制点）。 |
Ele |
等几何单元 |
GauPts |
高斯积分点的位置和权重。 |
dRu |
NURBS基函数对第一个参数方向（u）的导数。 |
dRv |
NURBS基函数对第二个参数方向（v）的导数。 |
参数空间->物理空间的雅可比矩阵推导如下：
根据定义
\[
\mathbf{J} =
\frac{\partial (x, y)}{\partial (\xi, \eta)} =
\begin{bmatrix}
\frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\
\frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta}
\end{bmatrix}
\]
上式的x是高斯点的x，不是控制点，根据下式可知，高斯点的坐标是由控制点插值得到
\[
\mathbf{x} =
\begin{bmatrix}
x \\
y
\end{bmatrix}
=
\sum_{i=1}^{n} N_i(\xi, \eta)
\begin{bmatrix}
x_i \\
y_i
\end{bmatrix}
\]
故每个元素对参数坐标的导数如下
\[
\frac{\partial x}{\partial \xi} = \sum_{i=1}^{n} \frac{\partial N_i}{\partial \xi} \cdot x_i
\]
\[
\frac{\partial x}{\partial \eta} = \sum_{i=1}^{n} \frac{\partial N_i}{\partial \eta} \cdot x_i
\]
\[
\frac{\partial y}{\partial \xi} = \sum_{i=1}^{n} \frac{\partial N_i}{\partial \xi} \cdot y_i
\]
\[
\frac{\partial y}{\partial \eta} = \sum_{i=1}^{n} \frac{\partial N_i}{\partial \eta} \cdot y_i
\]
带入J可知
\[
\mathbf{J} =
\begin{bmatrix}
\frac{\partial N_1}{\partial \xi} & \cdots & \frac{\partial N_n}{\partial \xi} \\
\frac{\partial N_1}{\partial \eta} & \cdots & \frac{\partial N_n}{\partial \eta}
\end{bmatrix}
\cdot
\begin{bmatrix}
x_1 & y_1 \\
x_2 & y_2 \\
\vdots & \vdots \\
x_n & y_n
\end{bmatrix}
\]
体现在代码中只有两行
1 2 | |
接下来为了计算位移-应变矩阵B，需要基函数对物理坐标的导数，使用链式法则如下：
\[
\frac{\partial N}{\partial x} =
\frac{\partial N}{\partial \xi} \cdot \frac{\partial \xi}{\partial x} +
\frac{\partial N}{\partial \eta} \cdot \frac{\partial \eta}{\partial x}
\]
\[
\frac{\partial N}{\partial y} =
\frac{\partial N}{\partial \xi} \cdot \frac{\partial \xi}{\partial y} +
\frac{\partial N}{\partial \eta} \cdot \frac{\partial \eta}{\partial y}
\]
转变为矩阵形式：
\[
\begin{bmatrix}
\frac{\partial N}{\partial x} \\
\frac{\partial N}{\partial y}
\end{bmatrix}
=
\begin{bmatrix}
\frac{\partial \xi}{\partial x} & \frac{\partial \eta}{\partial x} \\
\frac{\partial \xi}{\partial y} & \frac{\partial \eta}{\partial y}
\end{bmatrix}
\cdot
\begin{bmatrix}
\frac{\partial N}{\partial \xi} \\
\frac{\partial N}{\partial \eta}
\end{bmatrix}
\]
可以发现
\[
\mathbf{J}_1 =
\frac{\partial (x, y)}{\partial (\xi, \eta)}
\implies
\mathbf{J}_1^{-1} =
\frac{\partial (\xi, \eta)}{\partial (x, y)}
\]
故
\[
\begin{bmatrix}
\frac{\partial N}{\partial x} \\
\frac{\partial N}{\partial y}
\end{bmatrix}
=
\mathbf{J}_1^{-1} \cdot
\begin{bmatrix}
\frac{\partial N}{\partial \xi} \\
\frac{\partial N}{\partial \eta}
\end{bmatrix}
\]
在代码中体现如下
1 | |
从母空间到参数空间的转换如下：
\[
\begin{bmatrix}
u \\
v
\end{bmatrix}
=
\begin{bmatrix}
\frac{u_{\text{max}} - u_{\text{min}}}{2} & 0 \\
0 & \frac{v_{\text{max}} - v_{\text{min}}}{2}
\end{bmatrix}
\begin{bmatrix}
\xi \\
\eta
\end{bmatrix}
+
\begin{bmatrix}
\frac{u_{\text{max}} + u_{\text{min}}}{2} \\
\frac{v_{\text{max}} + v_{\text{min}}}{2}
\end{bmatrix}
\]
故雅可比矩阵\(J_2\)为
\[
J_2 =
\frac{\partial(u, v)}{\partial(\xi, \eta)} =
\begin{bmatrix}
\frac{u_{\text{max}} - u_{\text{min}}}{2} & 0 \\
0 & \frac{v_{\text{max}} - v_{\text{min}}}{2}
\end{bmatrix}
\]
整体代码注释如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 | |
Stiff_Ass2D
函数接口如下：
组装所有的局部刚度矩阵
1 | |
| 变量名 | 含义说明 |
|---|---|
KE |
每个单元的局部刚度矩阵（cell数组） |
CtrPts |
控制点结构体，含有总数 CtrPts.Num |
Ele |
单元结构体，包含每个单元关联的控制点编号等信息 |
Dim |
空间维度（通常为2，表示二维问题） |
Dofs_Num |
全局自由度总数（= 控制点个数 × 维度） |
整体代码注释如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 | |
Warning
全局刚度矩阵的规模为 \(自由度 \times 控制点个数 \times自由度 |times 控制点个数\)，sparse会自动将相同行列的下标处的值相加，刚度矩阵表示的是控制点自由度间的传递关系
solving
函数接口如下：
求解器，根据 \(U=K/F\)求解全局位移
1 | |
| 参数名 | 说明 |
|---|---|
CtrPts |
控制点信息，包括控制点总数 CtrPts.Num 等 |
DBoudary |
边界条件信息，包含被固定的控制点编号 |
Dofs |
自由度信息，含总自由度数 Dofs.Num、固定自由度等 |
K |
全局刚度矩阵，大小为 2N × 2N，其中 N 为控制点数 |
F |
全局载荷向量，大小为 2N × 1 |
BoundCon |
边界条件编号，表示具体结构问题类型 |
| 输出参数 | |
U |
全局自由度的位移解，大小为 2N × 1 |
整体代码注释如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 | |
3.7 sensitivity analysis
本章进行灵敏度分析
下式为目标函数对DDF的导数
\[
\begin{cases}
\frac{\partial J}{\partial \mathcal{X}} = -\frac{1}{2} \int_{\Omega} \varepsilon(\mathbf{u})^T \left( \gamma \mathcal{X}^{\gamma - 1} \mathbf{D}_0 \right) \varepsilon(\mathbf{u}) \, d\Omega \\
\frac{\partial G}{\partial \mathcal{X}} = \frac{1}{|\Omega|} \int_{\Omega} v_0 \, d\Omega
\end{cases}
\]
由于DDF是由平滑后的密度值插值得到的（由4.5章可知），
\[
\mathcal{X}(\xi,\eta) = \sum_{i,j} R^{p,q}_{i,j}(\xi,\eta) \cdot \tilde{\rho}_{i,j}
\]
故通过链式法则DDF对原密度值通过有如下表达式：
\[
\frac{\partial \mathcal{X}}{\partial \rho_{i,j}} = \frac{\partial \mathcal{X}}{\partial \tilde{\rho}_{i,j}} \cdot \frac{\partial \tilde{\rho}_{i,j}}{\partial \rho_{i,j}} = R^{p,q}_{i,j}(\xi,\eta) \cdot \psi(\rho_{i,j})
\]
则最终目标函数对设计变量的表达式如下：
\[
\begin{cases}
\frac{\partial J}{\partial \rho_{i,j}} = -\frac{1}{2} \int_{\Omega} \varepsilon(\mathbf{u})^T \gamma \mathcal{X}^{\gamma - 1} R^{p,q}_{i,j}(\xi,\eta) \psi(\rho_{i,j}) \mathbf{D}_0 \varepsilon(\mathbf{u}) \, d\Omega \\
\frac{\partial G}{\partial \rho_{i,j}} = \frac{1}{|\Omega|} \int_{\Omega} R^{p,q}_{i,j}(\xi,\eta) \psi(\rho_{i,j}) v_0 \, d\Omega
\end{cases}
\]
整体代码注释如下：
对于体积对设计变量的导数
\[
\frac{\partial G}{\partial \rho_{i,j}} = \frac{1}{|\Omega|} \int_{\Omega} v_0 \, d\Omega = \frac{1}{|\Omega|} \int_{\Omega} \frac{\partial \mathcal{X}(\xi, \eta)}{\partial \rho_{i,j}} \, d\Omega
\]
而由上面推导可知：
\[
\frac{\partial \mathcal{X}(\xi, \eta)}{\partial \rho_{i,j}} = R^{p,q}_{i,j}(\xi, \eta) \cdot \psi(\rho_{i,j})
\]
故
\[
\frac{\partial G}{\partial \rho_{i,j}} = \frac{1}{|\Omega|} \int_{\Omega} R^{p,q}_{i,j}(\xi, \eta) \cdot \psi(\rho_{i,j}) \, d\Omega
\]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | |
3.8 OC法
函数接口：
1 | |
| 参数名 | 含义 |
|---|---|
X.CtrPts |
控制点处的当前密度（设计变量），为列向量 |
R |
NURBS 基函数矩阵，用于将控制点密度映射到高斯点密度 |
Vmax |
最大允许体积分数（约束材料使用量） |
Sh |
Shepard 平滑矩阵，用于平滑设计变量 |
Hs |
Shepard 平滑归一化因子，避免除以0 |
dJ_dp |
目标函数对控制点密度的导数（柔度灵敏度） |
dv_dp |
体积约束函数对控制点密度的导数（体积灵敏度） |
函数整体注释如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | |
3.9 Plot_Data and Plot_Topy
Plot_Data：生成用于绘图的图像窗口布局参数（Pos），计算用于绘图的密度场坐标网格（DenFied）
1 | |
| 参数名 | 类型 | 含义说明 |
|---|---|---|
Num |
向量 | NURBS 网格在两个方向的单元数量，如 [20, 5] 表示 20×5 网格 |
NURBS |
结构体 | NURBS 几何结构，包含控制点、阶数、节点矢量等 |
DenFied |
结构体 | 绘图所需的密度场数据，包括基函数、物理坐标网格等 |
Pos |
结构体 | 图像窗口位置参数，用于多子图布局（p1 到 p4） |
Plot_Topy：等几何拓扑优化结果的可视化函数，绘制了设计变量分布、密度场和结构拓扑图等图形。
1 | |
| 参数名 | 类型 | 说明 |
|---|---|---|
X |
结构体 | 包含控制点密度 CtrPts、高斯点密度 GauPts、输出 DDF 等 |
GauPts |
结构体 | 高斯点的物理坐标 PCor 等 |
CtrPts |
结构体 | 控制点坐标 Cordis，数量等 |
DenFied |
结构体 | 绘图密度场所需信息，包括 Ux, Uy, N, id |
L |
数值 | 结构在 X 方向的长度（用于设定绘图坐标轴范围） |
W |
数值 | 结构在 Y 方向的宽度 |
Pos |
结构体 | 图形窗口的布局位置（由 Plot_Data 函数生成） |
四、实验结果
本章展示五种不同的模型，每种模型分别展示五副图，依次为：
- 控制点处的密度
- 高斯点处的密度
- 密度分布函数 DDF
- 高斯点处密度大于 0.5 的二维投影
- DDF 的等值线，近似表示结构拓扑边界
4.1 悬臂梁
4.2 MBB梁
4.3 Michell 型结构
4.4 L 型梁
4.5 四分之一圆环
