简介这是一份基于MATLAB的等几何分析IGA程序包nliga_v1.3面向计算力学与数值仿真初学者旨在解决传统有限元因网格近似导致的几何精度损失问题。资源共240个文件以194个m脚本为主辅以cc/c/h底层源码、msh网格文件、mexw64预编译组件及说明文档压缩包约12.34MB目录结构清晰便于按模块研读。程序涵盖NURBS曲线曲面基础函数库、IGA模型构建、线性与非线性ODE/PDE求解器、位移应力后处理可视化并配有分步示例与教程。已有1411人学习下载适合希望掌握等几何分析原理、NURBS数学概念以及MATLAB数值实现的科研与工程人员借助完整代码与案例可快速上手实践。 做这个 matlab 编制的 IGA等几何分析法程序起因其实很朴素我当时在做一个结构分析的小项目用传统有限元算完以后总觉得不对劲几何模型和计算网格之间隔着一层翻译而 IGA 恰好能把这层翻译省掉。等几何分析这个词听起来高级但说白了就是把 CAD 里描述几何的那套数学语言NURBS 样条直接拿来做数值分析的基函数网格就是几何几何就是网格。这套 matlab 程序就是我当时从零搭起来的一个可运行、可扩展的 IGA 小框架适合正在学等几何分析、或者想从有限元切换思路的朋友参考。这个程序能做的事情很直接给定一组 NURBS 控制点和节点向量它可以自动生成基函数计算单元刚度矩阵组装成全局方程施加边界条件最后求解并可视化位移场。你不需要外部依赖只要装了 matlab 就能跑而且因为矩阵运算全部用 matlab 原生语法代码量比 C 版本小很多特别适合用来验证算法思路。1. IGA 的思路与 matlab 选型考量1.1 从有限元到等几何分析到底改了什么传统有限元分析的第一步是网格划分把几何体离散成一个个小单元再用分片多项式去逼近真实形状。问题在于CAD 里描述的圆弧、样条曲面一旦离散成三角形或四边形就产生了几何误差。细化网格可以减小误差但永远是在逼近不是精确表达。IGA 的思路是既然 CAD 里用 NURBS 描述曲线曲面而有限元又需要基函数那为什么不直接用 NURBS 基函数作为形函数这样几何模型和计算模型完全统一粗网格时就是精确的几何形状细化时也不再需要重新建立几何只需要做节点插入。这个理念最早是 Hughes 团队在 2005 年提出的到现在已经成了数值方法研究里的一个成熟方向。用 IGA 还有一个容易忽略的好处因为 NURBS 基函数是高阶连续的比传统有限元的 C0 连续更光滑所以在涉及高阶导数的问题比如 Kirchhoff 板壳、梯度弹性里IGA 有天然优势。我的这个程序最初就是为了算实体结构但也顺手验证过梁、板问题效果都不错。1.2 为什么偏偏用 matlab 来写而不是 C 或 Python选 matlab 有几个非常现实的理由。第一矩阵操作太方便。IGA 的核心计算是基函数求值、雅可比映射、刚度矩阵积分这些在 matlab 里就是几行矩阵乘法和reshape的事情。第二调试可视化一体。plot、surf、patch可以直接把控制网格、基函数、位移场画出来这对理解 IGA 太重要了。你光看公式很难想象一个二维 NURBS 基函数长什么样但画出来立刻明白。第三也是比较实际的一点我的很多前期数据都是 mat 文件用 matlab 处理起来无缝衔接。后来虽然有人建议我转 Python但我发现 IGA 在 Python 里写虽然语法更现代但科学计算库的组合反而没有 matlab 紧凑。对个人研究者来说matlab 的所见即所得省掉了很多工程化精力。当然matlab 的短板是循环效率低。IGA 的组装过程往往涉及单元循环如果直接写逐点双重循环计算规模一大就会卡死。我的解决办法是用向量化 预分配来优化后面会详细讲。2. 程序整体设计与模块拆解2.1 模块划分一个 IGA 程序应该有哪些部分我按照功能把程序分成了六个模块几何定义、基函数计算、单元映射、刚度组装、边界条件、后处理可视化。每个模块单独成函数这样后续增加维度、增加细化算法都方便。几何定义模块读入 NURBS 曲线/曲面的控制点、权重、节点向量。基函数计算模块给定节点向量、阶数、参数点计算 NURBS 基函数及其导数。单元映射模块把参数空间的高斯点映射到物理空间计算雅可比矩阵。刚度组装模块计算局部刚度矩阵并组装到全局矩阵。边界条件模块处理 Dirichlet 边界固定约束因为是直接基于自由度施加。后处理模块把求解结果与控制网格叠加显示。这些模块里最核心的就是基函数计算它决定了一切。NURBS 基函数是 B 样条基函数的加权归一化而 B 样条基函数又有递推定义Cox-de Boor 递推。在 matlab 里递推可以用循环实现也可以写成递归函数但为了效率我建议用循环并把一维的情况做成一个独立的bspline_basis函数。2.2 数据结构控制点、节点向量与权重对于一个二维 IGA 程序需要两个节点向量( U ) 和 ( V )分别对应参数坐标 ( \xi ) 和 ( \eta )控制点是一个 ( n \times m \times 2 ) 的三维数组权重是对应的二维数组。我当时设计了这样的数据结构U [0, 0, 0, 0.5, 1, 1, 1]; % 节点向量阶数 p2有三个重复端节点 P zeros(n, m, 2); % 控制点坐标 w ones(n, m); % 权重 p 2; q 2; % 两个方向的样条阶数节点向量是整个程序最容易出错的地方。IGA 里节点向量不一定要均匀而且 open 节点向量首尾节点重复 ( p1 ) 次是标准做法因为它保证基函数在端点插值这样施加边界条件时可以直接对应到控制点自由度。如果你用均匀节点向量但没有重复端点则首尾边界上无法像有限元那样简单固定需要额外处理新手很容易在这里踩坑。2.3 与有限元程序的差异点一开始我是照搬有限元的矩阵组装思路但后来发现两者有一点必须调整有限元里单元是显式网格每个单元有自己的节点编号而 IGA 里单元是节点向量区间 ([U_i, U_{i1}])对应一组非零基函数的编号。所以 IGA 的单元-自由度连接关系是动态计算的不能像有限元那样预先生成单元节点矩阵。我的做法是写一个iga_dof_mapping函数遍历每个非空单元区间找到该区间上非零基函数的全局编号存成 cell 数组。组装时直接按每个单元的dofs列表来定位全局矩阵位置。这样程序就非常灵活换阶数、换细化程度都不用改主循环。% 伪代码查找非零基函数编号 function [indices, spans] find_nonzero_basis(U, p, nev) % nev: 需要求值的参数点数组 % 返回每个点的非零基函数起始索引 spans zeros(size(nev)); for i 1:length(nev) spans(i) find_span(U, p, nev(i)); % 二分查找 end indices spans - p : spans; % 非零基函数编号 end实际代码里find_span是基础工具负责在节点向量中定位参数点属于哪个节点区间。这个函数网上有很多现成实现但要小心边界情况当参数点恰好等于节点向量最后一个值时要特殊处理否则索引会越界。3. 核心环节实现与关键逻辑3.1 NURBS 基函数及其导数的递推计算如果你已经熟悉 B 样条基函数的 Cox-de Boor 递推那 NURBS 就只是加了一个权重归一化计算顺序是用递推公式计算 B 样条基函数值 ( N_{i,p} ) 及其一阶导/二阶导。计算权重和 ( W(\xi) \sum N_{i,p} w_i )。计算 NURBS 基函数 ( R_i N_{i,p} w_i / W(\xi) )。对于导数用除法求导法则 [ R \frac{N w W - N w W}{W^2} ]我写的函数支持输出 0 阶、1 阶、2 阶导因为后续算刚度矩阵时要用到一阶导算曲率问题时需要二阶导。由于是用循环实现的建议提前把基函数值存入临时数组避免重复计算。效率上二维问题的基函数是两个一维基函数的张量积所以只需要分别算 ( R_i(\xi) ) 和 ( R_j(\eta) )再相乘即可这样大大减少了计算量。3.2 参数空间到物理空间的映射等几何分析和有限元最大的不同就在这里有限元的形函数是定义在参考单元上的而 IGA 的形函数直接定义在几何参数空间里。因此计算单元刚度矩阵时需要把参数空间的高斯积分点映射到物理空间。对于二维问题这个映射是[ \mathbf{x}(\xi,\eta) \sum_{i}\sum_{j} R_{i,j}(\xi,\eta) \mathbf{P}_{i,j} ]映射的雅可比矩阵为[ \mathbf{J} \begin{bmatrix} x_\xi y_\xi \ x_\eta y_\eta \end{bmatrix} ]在程序里我先用基函数对 ( \xi ) 和 ( \eta ) 的偏导数分别得到 ( x_\xi, x_\eta, y_\xi, y_\eta )然后组合成 2x2 矩阵。如果处理三维实体曲面雅可比就是 3x2 矩阵刚度积分时需要用 ( \det(\mathbf{J}^T \mathbf{J}) ) 代替 ( \det(\mathbf{J}) )。我当时在二维程序里统一用一个detJ变量三维版本里换成了sqrt(det(J*J))。3.3 高斯积分点的处理IGA 单元的高斯积分和有限元完全相同但有一点须注意每个非空节点区间内基函数的连续性可能不同因此需要足够多的积分点。一般来说每方向取阶数 ( p1 ) 个高斯点足够精确。但如果你把基函数细化以后单元尺寸变小积分点不需要增加因为基函数在每个单元内的多项式性质不变。我在程序中定义了gauss_quadrature函数返回高斯点和权重然后对所有单元循环。为了提高效率我预先计算了所有单元在高斯点的基函数值而不是在每个单元里重复调用递推函数。这个优化让组装时间降低了大约 60%在循环密集型 matlab 代码里效果非常明显。3.4 刚度矩阵组装与自由度处理以二维线弹性问题为例单元刚度矩阵的通用形式是[ \mathbf{K}^e \int_{\Omega^e} \mathbf{B}^T \mathbf{D} \mathbf{B} , \det(\mathbf{J}) , d\xi d\eta ]其中 ( \mathbf{B} ) 是应变-位移矩阵由基函数导数组成( \mathbf{D} ) 是材料本构矩阵。组装过程与有限元一样通过单元自由度映射dofs将局部矩阵累加到全局矩阵。边界条件的施加我采用最直接的方法先求解再把固定自由度上的位移置零并修改对应的行列。如果固定位移不为零就需要处理右端项的修正。更稳妥的做法是划零置一法但它在 matlab 稀疏矩阵上反而比较慢所以我用了一个简便做法% 施加 Dirichlet 边界条件固定位移 free_dofs setdiff(1:ndof, fixed_dofs); F_modified F(free_dofs) - K(free_dofs, fixed_dofs) * u_fixed; u zeros(ndof, 1); u(fixed_dofs) u_fixed; u(free_dofs) K(free_dofs, free_dofs) \ F_modified;这里最关键的是搞清楚自由度编号与控制点编号的对应关系。对于二阶张量问题每个控制点有 ( d ) 个自由度2D 问题里 ( d2 )所以全局自由度编号是dof (control_point_index - 1)*d component。这块映射错了程序会报奇异性或者结果完全错乱。我的经验是先写一个输出dofs和固定点的联动函数然后单独打印前 20 个映射关系肉眼检查一遍再继续。3.5 后处理直接绘制 NURBS 几何和位移场后处理是 IGA 比传统 FEM 漂亮的地方。因为粗网格的控制点并不在物理曲线上不能简单画控制点连线来代表变形。正确的做法是在参数空间里采样大量点用 NURBS 基函数映射到物理空间。我写了evaluate_nurbs_surface函数对给定的参数网格逐一求值然后与位移场叠加% 绘制变形后的几何 X_deformed X reshape(U_x, size(X)); surf(X_deformed, Y_deformed, 0*X_deformed)如果想把控制网格也画出来直接连控制点即可。控制网格和真实几何之间的关系可以非常直观地展示等几何的含义——控制网格像是一个骨架而几何是骨架撑起来的光滑表面。4. 实操过程与参数调试4.1 从二维悬臂梁开始跑通流程我建议第一次跑程序不要用复杂模型用一个单位方形板或者悬臂梁就够了因为解析解容易对比。我用的例子是 2x1 矩形板左端固定右端施加均布拉力材料采用单位杨氏模量泊松比 0.3平面应力。这种情况下基函数阶数 ( pq2 )控制点分别是 4x3 个点。节点向量设为U [0 0 0 1 1 1]; V [0 0 0 1 1 1];这样只有一个单元但因为是二次 NURBS它可以精确表示线性的位移场所以解的结果和理论值几乎完全一致误差在 1e-14 量级。跑通这个例子后再尝试细化在 ( U ) 方向插入节点 0.5变成U [0 0 0 0.5 1 1 1]程序会自动把单元一分为二不需要重新建立几何。这就是 IGA 的一个关键优势细化是插入节点而不是重新网格化。4.2 如何选择阶数和节点向量分布很多初学者一上来就选三阶、四阶觉得阶数越高越精确。但实际上阶数越高非零基函数的跨度越大组装成本越高而且边界附近容易出现振荡。我的实际经验是二维静力分析二阶已经足够如果涉及光滑性或高阶导数再用三阶。节点向量分布取决于几何有圆角或曲率变化大的区域节点要加密但不是均匀加密而是按照曲率来插。插入节点也有讲究不能随便在某个位置插入IGA 的标准做法是等分参数区间或按曲率自适应。我当时写过一个knot_insert函数它按照 Hughes 书里的算法实现核心是新旧控制点的线性组合。这个算法一两行说不清但网上有现成 matlab 代码可参考抄过来自己调试一遍能加深理解。4.3 收敛性验证与误差分析我跑完后最喜欢做的事是做一些收敛性对比。对于二维弹性问题理论上的收敛率是 ( O(h^{p1}) )即阶数越高、网格越细误差下降越快。我用 ( L_2 ) 误差范数来评价L2err sqrt(sum(((u_num - u_exact).^2) .* area_per_point));这里要注意采样点要足够密否则误差算不准。我通常在一个单元里取 20x20 个采样点虽然效率略低但结果可信。对比下来二阶 IGA 和一阶线性有限元在相同自由度下的精度差异非常显著——IGA 往往用更少的自由度就能达到相同精度这就是高阶连续带来的红利。5. 常见问题与排查技巧实录5.1 刚度矩阵奇异或求解失败最常见的原因是自由度映射错误或约束不足。如果刚度矩阵零特征值数量明显多于刚体模态数基本就是固定边界没施加对。我的建议是把刚体模态导出来检查对二维实体问题不加约束时应该有 3 个零特征值两个平动、一个转动。如果多了就说明个别自由度没有被连接或约束问题多半出在setdiff之后u的维度对不上。另一个容易忽略的原因是重复的控制点或权重为零。NURBS 中如果某个控制点权重为 0基函数为 0那这个自由度相当于不存在但你又给它施加约束就可能造成矩阵不满秩。我在程序里加了检查any(w eps)直接报错避免这种隐性 bug。5.2 基函数在节点处数值跳变如果参数点恰好落在重复节点上基函数可能存在不连续导数。你的积分点要避开这些位置所以我在高斯积分时使用gauss_rule的内部点而不是单元端点。但如果需要精确评估节点处的值例如固定位移点建议取节点向量区间内部的一点来逼近或者在重复节点处直接采用极限值。5.3 程序运行很慢怎么优化matlab 循环慢IGA 又天然有多个嵌套循环单元循环、高斯点循环、基函数循环。我的优化策略有三条预计算基函数数值把二维基函数张量积事先生成好避免在每个单元里重复调用递推。用稀疏矩阵存储全局刚度矩阵直接用sparse(i, j, s, m, n)从行列索引向量构造而不是在一个个循环里赋值。把最内层的小循环替换为矩阵乘法。例如计算局部刚度矩阵时将 B 矩阵写成块矩阵然后一次性乘上 D 矩阵而不是逐个分量乘。这样优化后一万自由度左右的算例在普通笔记本上几秒钟就能完成性能已经可以接受。5.4 边界条件施加的自由度编号错位边界条件不值得炫技最容易错的是编号。IGA 的控制点编号、自由度编号和你在几何定义时的排序不同就全乱了。我的血泪教训是写一个show_index调试函数把控制点坐标和自由度编号一起打印出来同时可视化控制点编号与几何对应一眼就能发现哪个编号错位。比如二维问题控制点按i、j两个方向索引自由度是 ( (i,j) ) 处有 ( u_x, u_y ) 两个。固定左端时要先把i1的索引取出来再映射到自由度。如果你用的是sub2ind来生成自由度要记得维度顺序是[dim, ny, nx]否则就会把所有点都映射错。最后再分享一个小技巧在写 IGA 程序时不要一次性追求功能完备先从一个单元、一个方向开始验证基函数求值和刚度矩阵组装正确再逐步扩展到二维、三维。我当时的调试次序是先画 NURBS 几何再输出基函数值对照书上的数值逐一验证然后算悬臂梁位移最后才做细化。这个次序虽然慢但每一步出错都能快速定位。实际上IGA 程序真正难的地方不在于某个函数有多复杂而在于所有环节之间的一致性——从节点向量到基函数再到自由度映射任何一处不匹配都会导致整个求解崩盘。希望这套 matlab 程序能帮你少走一些弯路如果你在跑的过程中遇到什么问题欢迎对照上面的排查思路逐项检查。本文还有配套的精品资源点击获取