简介这份资源是面向初学者的MATLAB等几何分析IGA程序包nliga_v1.3用于在MATLAB环境中学习并实践将CAD参数化几何直接用于数值模拟的方法。包内共240个文件以194个.m脚本与函数文件为主体辅以12个.cc与7个.c底层NURBS基函数实现、11个.msh网格数据、6个.mexw64编译接口及少量头文件与说明文档压缩包约12.34MB结构上覆盖基础函数库、模型构建、求解器与后处理等模块。已有1415人学习下载说明其在入门群体中具有一定参考价值。读者可借助其中的NURBS曲线曲面构造、B样条基函数计算与示例案例理解IGA基本原理并动手运行简单分析逐步掌握MATLAB编程与结构、流体等工程问题数值模拟的衔接方式。1. 拿到一份 MATLAB 版 IGA 程序先别急着跑前阵子帮一个做结构优化的朋友看代码他丢过来一个压缩包说“这是别人用 MATLAB 写的等几何分析程序我跑不通”。我打开一看目录里全是.m文件没有 README没有示例脚本函数名还都是buildK、assemble、refine这种通用词。这种“裸源码”在工程圈太常见了——作者自己跑得通换台机器就玄学报错。等几何分析IGA本身不复杂用 NURBS 基函数代替传统有限元的拉格朗日插值让几何模型和求解空间共用一套样条表示省掉网格划分这一步。这份 MATLAB 程序就是把 IGA 的完整链路——控制点、节点向量、基函数求值、刚度矩阵组装、边界条件施加、线性方程组求解——用脚本串起来。适合两类人学过有限元但没碰过样条基的学生以及想把 IGA 嵌进自己优化流程的工程师。下面按“能干什么、怎么跑、怎么改”拆开讲。2. 程序骨架拆解从节点向量到刚度矩阵的完整链路2.1 目录与函数分工拿到源码先别逐行读先看文件命名和调用关系。这类 IGA 程序通常分四层几何层控制点坐标、权重、节点向量、基函数层B 样条/NURBS 基函数及其导数、组装层单元循环、高斯积分、刚度矩阵、求解层边界条件、线性求解、后处理。我见过的典型文件清单长这样文件职责调用方main.m参数配置、流程串联用户直接运行knotVector.m生成/细化节点向量main.mbasisFuns.m单点基函数及导数求值组装层buildGeometry.m控制点、权重、连接关系main.massembleStiffness.m单元循环 高斯积分main.mapplyBC.m本质边界条件处理main.msolveSystem.m线性方程组求解main.m这个分层不是作者拍脑袋定的是 IGA 的数学结构决定的基函数只认节点向量和参数坐标组装只认基函数值和雅可比求解只认全局矩阵。你把basisFuns.m单独拎出来测输入一个参数坐标和节点向量输出基函数值对不上就是这里的问题不用往下查。常见做法是每个函数开头写一段注释说明输入输出维度但这份源码大概率没有——我一般会先补一个test_basis.m做单元测试拿一个已知的二次 B 样条节点向量[0 0 0 1 2 2 2]验证。2.2 节点向量与基函数IGA 的地基IGA 和传统 FEM 最大的区别就在这一层。传统 FEM 的形函数是固定的多项式IGA 的基函数由节点向量决定改变节点向量就改变了基函数的连续性和支撑域。看这段典型实现function [N, dN] basisFuns(i, u, p, U) % 计算第 i 个 p 次基函数在参数 u 处的值及一阶导数 % i: 基函数索引1-based % u: 参数坐标 % p: 次数 % U: 节点向量 N zeros(p1, 1); dN zeros(p1, 1); left zeros(p1, 1); right zeros(p1, 1); N(1) 1.0; for j 1:p left(j1) u - U(i1-j); right(j1) U(ij) - u; saved 0.0; for r 0:j-1 temp N(r1) / (right(r2) left(j-r1)); N(r1) saved right(r2) * temp; saved left(j-r1) * temp; end N(j1) saved; end % 导数计算略按同样递推结构展开 end这段是 Cox-de Boor 递推的原地版本left和right存的是参数点到左右节点的距离saved用来滚动保存中间量。参数说明里最容易翻车的是索引MATLAB 从 1 开始而文献里的公式从 0 开始U(i1-j)这种写法就是在做偏移。如果你改次数p节点向量长度必须满足length(U) n p 1n 是控制点数否则递推直接越界。我一般会在main.m里加一句断言assert(length(U) size(P,1) p 1, 节点向量长度与控制点数不匹配);这句断言能挡掉八成“跑不通”的问题——很多人改了几何参数忘了同步节点向量。2.3 组装层单元循环与高斯积分组装层的逻辑和标准 FEM 几乎一样区别只在形函数来源。核心循环长这样function K assembleStiffness(P, W, U, V, p, q, ngp) % P: 控制点坐标 (n x 2) % W: 权重 (n x 1) % U, V: 两个参数方向的节点向量 % p, q: 两个方向的次数 % ngp: 每个方向的高斯点数 n size(P, 1); K sparse(n, n); [gp, gw] gaussQuad(ngp); % 高斯点与权重 % 遍历参数域单元简化单单元情形 for a 1:ngp for b 1:ngp u 0.5*(1gp(a)); % 映射到 [0,1] v 0.5*(1gp(b)); [R, dRdu, dRdv] nurbsBasis(u, v, U, V, p, q, W); J dRdu * P; % 雅可比矩阵 detJ det(J); dRdx J \ [dRdu dRdv]; % 物理域导数 K K (dRdx * dRdx) * detJ * gw(a) * gw(b); end end end几个关键点sparse预分配是必须的控制点数上千后稠密矩阵直接爆内存J \ [dRdu dRdv]是解雅可比线性系统求物理域导数比显式求逆稳detJ必须为正出现负值说明单元参数方向反了或者控制点顺序错了。高斯点数ngp的经验值是ceil((pq2)/2)二次基函数取 3 点足够取多了浪费取少了积分不精确——这个坑我在板壳单元上踩过刚度矩阵不对称查了半天才发现是积分欠采样。2.4 边界条件与求解本质边界条件的特殊处理IGA 最反直觉的地方在这里NURBS 基函数不满足 Kronecker delta 性质控制点上的值不等于物理场在该点的值。所以传统 FEM 里“把对应自由度置零”的做法不能直接用。常见处理有三种罚函数法、Nitsche 法、控制点插值法。这份源码大概率用的是最粗暴的罚函数function [K, F] applyBC(K, F, fixedDofs, penalty) % fixedDofs: 被约束的自由度索引 % penalty: 罚因子通常取 1e8 ~ 1e12 for i 1:length(fixedDofs) d fixedDofs(i); K(d, d) K(d, d) penalty; F(d) F(d) penalty * 0; % 齐次边界 end end罚因子取太小约束不紧取太大条件数爆炸。我一般取1e10量级然后检查残差norm(K*x - F)/norm(F)超过1e-6就说明罚因子不合适或者基函数有问题。如果你要做非齐次边界比如给定位移F(d)那行要改成penalty * prescribedValue但注意这个值不是控制点值得先做插值换算——这是二次开发时最容易写错的一行。3. 跑通第一遍环境、参数与验证3.1 环境准备与最小运行路径MATLAB 版本不用太新R2016b 以上都行因为代码里大概率没用到字符串数组或隐式扩展这些新特性。不需要额外工具箱sparse、det、\都是基础功能。跑之前先确认三件事当前目录切到源码根目录否则函数找不到、main.m里的几何参数和节点向量自洽、输出目录有写权限。最小运行路径就是直接敲 main如果报“未定义函数或变量”八成是路径问题如果报“索引超出数组范围”去查节点向量长度断言如果跑完但位移全零查边界条件索引是不是 1-based 对上了。我习惯在main.m末尾加一段自检residual norm(K*x - F) / norm(F); fprintf(相对残差: %.3e\n, residual); assert(residual 1e-6, 求解残差过大检查罚因子或基函数);这段自检能帮你区分“跑完了”和“跑对了”——很多人看到不报错就以为成了结果位移场是错的。3.2 用已知解验证悬臂梁算例验证 IGA 程序最稳的办法是拿一个解析解已知的算例。二维悬臂梁端部受集中力Timoshenko 梁理论有解析解但 IGA 程序多半是平面应力或平面应变得用更简单的单向拉伸。左端固定右端均布载荷解析解是线性位移场。跑完提取右端控制点位移和理论值比% 理论值u F*L/(E*A) u_theory F_total * L / (E * thickness * height); u_num max(abs(x(1:2:end))); % 取 x 方向最大位移 err abs(u_num - u_theory) / u_theory; fprintf(相对误差: %.2f%%\n, err*100);误差在 1% 以内算正常超过 5% 就要查基函数次数够不够、高斯积分够不够、边界条件施加方式对不对。注意 NURBS 的权重要全设为 1 才退化成 B 样条如果源码里权重是随机给的先改成全 1 再验证。3.3 网格细化h-refinement 与 p-refinementIGA 的细化比 FEM 灵活因为可以只改节点向量不改几何。h-refinement 是插入节点p-refinement 是升阶k-refinement 是两者结合。源码里如果有refine.m大概率只实现了 h-refinement 的节点插入function [U_new, P_new, W_new] hRefine(U, P, W, p, numInsert) % 在参数域中点插入 numInsert 个节点 for k 1:numInsert u_new 0.5; % 简化固定中点 [U, P, W] knotInsert(U, P, W, p, u_new); end U_new U; P_new P; W_new W; end节点插入的公式Boehm 算法不难但索引偏移容易错。验证方法是插入前后几何不变在参数域取 100 个点用插入前后的基函数分别算物理坐标最大偏差应该小于1e-12。这个测试我每次都跑因为节点插入写错的话几何会“悄悄变形”后面所有结果都不可信。4. 避坑与排查五个真实翻车现场4.1 现象刚度矩阵不对称原因高斯积分点数不足或者基函数导数在单元边界上求值不连续。IGA 的基函数在单元内部是多项式高斯积分应该精确但如果ngp取少了积分误差会破坏对称性。解决把ngp提到pq2再检查norm(K - K)/norm(K)应该小于1e-14。4.2 现象位移场出现振荡原因罚因子过大导致条件数爆炸或者 NURBS 权重出现负值。解决先把罚因子降到1e8试再检查权重是否全正。NURBS 权重必须为正这是数学要求源码里如果从文件读权重读进来先assert(all(W 0))。4.3 现象细化后结果反而变差原因节点插入时控制点更新公式用错几何发生了微小变形。解决跑 3.3 节的几何不变性测试偏差超过1e-10就说明插入算法有问题。我见过一个版本把U(i1-j)写成了U(i-j)几何偏差 0.3%肉眼看不出来但应力集中区域误差翻倍。4.4 现象detJ出现负值原因控制点顺序反了或者参数方向映射搞反。解决检查P的排列是不是按参数方向递增二维问题里P通常按u快v慢排列。如果几何是从 CAD 导出的控制点顺序可能和程序约定不一致需要重排。4.5 现象换台机器结果不一样原因MATLAB 版本差异导致sparse求解器默认选项不同或者浮点累加顺序不同。解决在main.m开头固定随机种子如果有随机线性求解显式指定x K \ F而不是用迭代求解器。IGA 的刚度矩阵是对称正定的直接求解最稳。5. 二次开发把 IGA 嵌进自己的流程5.1 改造方向从单 patch 到多 patch这份源码大概率只支持单 patch一个 NURBS 面。实际工程里几何往往是多 patch 拼的改造的关键是 patch 间的连续性约束。常见做法是在公共边界上做控制点合并或者用 mortaring 方法。我一般先做最简单的公共边界上的控制点坐标强制相等用拉格朗日乘子施加。代码上就是在assembleStiffness之后加一个约束矩阵C解鞍点系统% 多 patch 连续性约束C*x 0 K_aug [K, C; C, zeros(size(C,1))]; F_aug [F; zeros(size(C,1),1)]; x_aug K_aug \ F_aug; x x_aug(1:size(K,1));这个改法的代价是系统不再正定得用\的通用求解。如果 patch 数量多鞍点系统会变大可以考虑静凝聚消去内部自由度。5.2 改造方向接外部几何数据从 CAD 导出 NURBS 数据通常是 IGES 格式MATLAB 没有原生读取函数。我一般让上游导出成简单的文本格式控制点坐标、权重、节点向量各一段用dlmread或readmatrix读。读进来第一件事是单位统一——CAD 常用毫米程序里如果按米算刚度差1e6倍结果全错。我习惯在读入后加一句scale 1e-3; % 毫米转米 P P * scale;然后跑 3.2 节的拉伸验证确认量纲对了再往下做。5.3 验证习惯每次改动都跑回归二次开发最怕“改一处坏三处”。我的习惯是维护一个小型回归集拉伸算例验量纲、悬臂梁算例验精度、细化算例验几何不变性。每次改完basisFuns或assembleStiffness三个算例全跑一遍残差和误差都打印出来对比。这套习惯是被坑出来的——有次改罚因子施加方式拉伸算例过了悬臂梁误差从 1% 跳到 8%查了一下午才发现是边界条件索引偏移了一位。从那以后我每次动边界条件相关代码都强制走一遍三个算例的回归不跳步。希望这份拆解能帮你少走几个弯路把这份 IGA 程序真正跑起来、改下去。本文还有配套的精品资源点击获取