简介：这份资源是面向初学者的MATLAB等几何分析（IGA）入门程序包nliga_v1.3，适合希望理解IGA原理、掌握NURBS数学概念并在MATLAB中实现数值模拟的学生与工程研究人员。压缩包共240个文件，约12.34MB，以194个.m脚本为主体，配合.c/.cc/.h底层源码、.mexw64编译模块、.msh网格数据及makefile等构建文件，覆盖NURBS基函数生成、模型构建、求解与后处理等环节。已有1415人学习下载。包内基础函数库可帮助读者理解控制点、权重与B样条曲面的构造方式，求解器与示例案例则便于对照运行、观察位移与应力结果，是衔接CAD几何与有限元分析的实用练习材料。

1. 一份能跑通的 MATLAB 等几何分析程序，到底解决了谁的痛点

如果你做过结构仿真，大概率经历过这种割裂：CAD 里是 NURBS 曲面，网格划分器里是线性四面体，两边几何永远对不上，光是把模型倒进求解器就得折腾半天，倒进去还得担心网格畸变。等几何分析（IGA）的思路很直接——既然 CAD 用 NURBS 描述几何，那就让 NURBS 基函数直接当形函数，几何和离散统一，省掉中间那层网格转换。这份 MATLAB 编制的 IGA 程序，核心价值就在这儿：它把 NURBS 基函数构造、控制点映射、刚度矩阵组装、边界条件施加这条链路用 MATLAB 脚本串起来了，适合做板壳、二维弹性体这类中小规模问题的教学复现和算法验证。MATLAB 的好处是矩阵运算天然贴合有限元组装，调试时能直接看基函数图像，比 C++ 改一行编译半天友好得多。适合谁？做结构仿真想入门 IGA 的研究生、需要快速验证新基函数或新积分方案的算法工程师，以及手头有 NURBS 几何但不想碰商业软件二次开发的人。

2. NURBS 基函数与等参映射：程序的地基怎么搭

2.1 从节点矢量到基函数值：程序里最不能省的一步

IGA 和传统有限元最大的差别，就在基函数不是插值型的拉格朗日多项式，而是 NURBS 基函数。程序里第一步永远是给定节点矢量（knot vector），比如一个 p 次、n 个控制点的一维情况，节点矢量长度是 n+p+1。MATLAB 里常见做法是先写一个openKnotVector函数生成 clamped 节点矢量，再写basisFuns用 Cox-de Boor 递推算基函数值。这段代码是整个程序的根，基函数算错，后面刚度矩阵全是废的。

function [N, dN] = basisFuns(i, u, p, U) % i: 当前span索引, u: 参数点, p: 次数, U: 节点矢量 % 返回 p+1 个非零基函数值 N 及其一阶导 dN N = zeros(1, p+1); dN = zeros(1, p+1); N(1) = 1; for j = 1:p left(j+1) = u - U(i+1-j); right(j+1) = U(i+j+1) - u; saved = 0; for r = 0:j-1 denom = right(r+2) + left(j-r+1); temp = N(r+1) / denom; N(r+1) = saved + right(r+2) * temp; saved = left(j-r+1) * temp; end N(j+1) = saved; end % 导数部分略，按标准递推补全 end

这段是 Cox-de Boor 的标准实现，left和right两个临时数组存的是当前参数点到左右节点的距离，denom是分母，saved用来滚动保存上一轮的值。参数说明：i必须由findSpan函数先算出来，不能随手给；U必须是单调不减的，重复度不能超过 p+1，否则分母为零直接 NaN。我一般会在调用前加一句assert(all(diff(U)>=0))做防御，血泪经验是节点矢量手写错一位，基函数图像看着还挺正常，但积分结果偏 3% 以上，查半天查不出来。

2.2 等参映射与雅可比：二维问题里最容易翻车的地方

一维基函数算完，二维就要做张量积。程序里通常用两个方向的节点矢量U、V和次数p、q，控制点排成(n+1)×(m+1)的网格。每个高斯积分点上，先分别算Nu、Nv，再N = kron(Nv, Nu)得到单元基函数向量。几何映射是x(ξ,η) = Σ N_A(ξ,η) * P_A，其中P_A是控制点坐标。雅可比矩阵J = [∂x/∂ξ ∂x/∂η; ∂y/∂ξ ∂y/∂η]，程序里必须用det(J)做积分权重，同时inv(J)把物理坐标导数转成参数坐标导数。

% 单元循环内，高斯点 (xi, eta) 处 [Ru, dRu] = basisFuns(spanU, xi, p, U); [Rv, dRv] = basisFuns(spanV, eta, q, V); R = kron(Rv, Ru); % 1 x (p+1)(q+1) dR = [kron(Rv, dRu); kron(dRv, Ru)]; % 2 x nloc J = dR * Pcoords; % 2 x 2, Pcoords 为单元控制点坐标 detJ = det(J); if detJ <= 0 error('雅可比非正，检查控制点顺序或节点矢量'); end dRdx = J \ dR; % 物理坐标导数

这里Pcoords是当前单元涉及的控制点坐标矩阵，行数等于局部基函数个数。detJ <= 0的判断是必须的，控制点顺序反了或者节点矢量重复度给错，雅可比会变负，积分出来的刚度矩阵不对称，求解器可能不报错但结果完全不对。常见做法是在单元循环里把detJ的最小值记下来，跑完打印一遍，如果接近零说明几何有退化，得回头查控制点。另一个坑是kron的顺序，kron(Rv, Ru)和kron(Ru, Rv)得到的基函数排列不一样，必须和后面组装时的自由度编号严格对应，否则刚度矩阵看着对称，实际是错位的。

2.3 刚度矩阵组装：稀疏存储和自由度映射

IGA 的组装逻辑和有限元几乎一样，差别在单元连接关系不是网格给的，而是由节点矢量的 span 推出来的。程序里一般先建一个elementConn矩阵，每行存一个单元的控制点全局编号，再建dofMap把控制点编号映射到自由度编号（二维弹性问题每个控制点两个自由度）。组装时用三元组(I, J, V)存稀疏矩阵，最后sparse(I,J,V,nDof,nDof)一次性生成。

nloc = (p+1)*(q+1); Ke = zeros(2*nloc, 2*nloc); for g = 1:nGauss % ... 算 R, dRdx, detJ, 权重 w ... B = zeros(3, 2*nloc); for a = 1:nloc B(:, 2*a-1:2*a) = [dRdx(1,a) 0; 0 dRdx(2,a); dRdx(2,a) dRdx(1,a)]; end Ke = Ke + B' * D * B * detJ * w; end % 全局组装 rows = dofMap(conn); cols = rows; I = [I; repmat(rows, 2*nloc, 1)]; ...

D是材料矩阵，平面应力问题就是那个含弹性模量和泊松比的 3×3 矩阵。B矩阵的构造顺序必须和dofMap一致，常见错误是dRdx第一行对应 ξ 方向导数、第二行对应 η 方向，但B里把 x/y 搞反，结果位移场方向整体旋转 90 度。组装完建议先检查norm(K - K')，IGA 刚度矩阵理论上对称，不对称说明组装索引有重复或遗漏。稀疏矩阵的nDof要按控制点总数乘自由度算，不是单元数，这一点和有限元不同，新手容易按单元节点数去开数组，直接越界。

3. 边界条件与求解：控制点自由度怎么约束才不出错

3.1 Dirichlet 边界条件的施加：强施加和罚函数的选择

IGA 里施加位移边界条件比有限元麻烦，因为 NURBS 基函数不是插值型的，控制点上的位移值不等于物理边界上的位移值。程序里常见两种做法：一是强施加，直接找到边界上对应的控制点，把对应自由度从方程组里消掉；二是罚函数法，在刚度矩阵对角上加一个大数。这份程序大概率用的是强施加，因为教学代码追求结果干净。

强施加的步骤：先根据参数域边界（比如 ξ=0 那条边）反算出对应的控制点索引，把这些控制点的自由度编号收集起来，然后对刚度矩阵做行列置零、对角置一、右端项置位移值。MATLAB 里可以用K(bcDof,:) = 0; K(:,bcDof) = 0; K(bcDof,bcDof) = speye(length(bcDof)); F(bcDof) = u_bc;。注意K是稀疏矩阵时，K(bcDof,:) = 0这种操作效率不高，控制点多了会明显变慢，常见做法是先转成sparse再操作，或者用逻辑索引一次性赋值。

提示：强施加时边界控制点的位移值不是随便给的，如果边界是曲线，控制点不一定落在物理边界上，直接给常数值会引入误差。教学代码里通常假设边界是直线或者控制点恰好落在边界上，实际工程问题要改用 Nitsche 方法或罚函数。

3.2 高斯积分阶数与求解器选择

IGA 单元内基函数是 p 次多项式，刚度矩阵被积函数次数是 2p，所以每个方向高斯点数至少取 p+1 才能精确积分。程序里常见做法是nGauss = p+1，但如果是变刚度或者几何非线性，得往上加。我一般会跑两组积分阶数对比，结果差超过 1e-6 就说明积分不够。求解器方面，MATLAB 反斜杠对稀疏对称正定矩阵会自动选 Cholesky，但 IGA 刚度矩阵条件数随次数升高恶化很快，p 大于 4 时建议显式用ichol加预条件，或者直接pcg带对角预条件。

% 求解并检查残差 K = (K + K')/2; % 强制对称，消除组装舍入误差 [L, U] = lu(K); % 小规模直接分解 u = U \ (L \ F); res = norm(K*u - F) / norm(F); fprintf('相对残差: %.3e\n', res);

(K+K')/2这一步不是多余，组装过程中浮点累加会让矩阵有 1e-15 量级的不对称，对称求解器可能报错。res超过 1e-10 就要警惕，要么是边界条件施加时行列置零破坏了对称性，要么是某个单元雅可比接近零导致矩阵病态。常见误用是直接u = K\F不检查残差，跑出来位移场看着合理，但应力云图在边界附近有异常振荡，回头查半天才发现是某个控制点自由度没约束住。

3.3 后处理：从控制点位移到物理场

解出来的是控制点位移，要画云图得在参数域上取密集采样点，用基函数插值回物理域。程序里一般写一个postProcess函数，在 ξ、η 方向各取 50 个点，算每个点的位移和应力，再用scatter或trisurf画出来。应力计算要注意，IGA 的应力在单元边界上连续，但程序里如果按单元循环算，边界点会被算两次，画云图时会出现一条缝。常见做法是采样点避开单元边界，或者用accumarray做平均。

nSample = 50; xiVec = linspace(0, 1, nSample); [Xi, Eta] = meshgrid(xiVec, xiVec); Ux = zeros(size(Xi)); Uy = zeros(size(Xi)); for a = 1:numel(Xi) % 找 span, 算 R, 插值位移 [~, spanU] = findSpan(...); R = ...; Ux(a) = R * uctrl(dofU); Uy(a) = R * uctrl(dofV); end trisurf(delaunay(Xi(:), Eta(:)), Ux(:));

采样密度nSample取 50 是经验值，p 高的时候要加到 80 以上，否则云图看着有棱角。delaunay在规则网格上效率不高，可以直接用surf加参数网格，但物理域不是矩形时得先映射。后处理这块最容易出的问题是位移插值用了错误的 span，参数点落在单元边界上时findSpan返回的索引可能和相邻单元不一致，导致云图在单元交界处跳变，解决办法是采样点稍微偏离边界，或者统一用左闭右开规则。

4. 避坑与排查：跑这份 IGA 程序时最常见的五类翻车

4.1 基函数图像正常但结果偏几个百分点

现象：跑一个悬臂梁算例，位移和理论解差 3% 到 5%，基函数画出来光滑漂亮，看不出问题。原因多半是节点矢量重复度给错，clamped 节点矢量首尾重复度应该是 p+1，中间节点重复度不超过 p，如果首尾只重复了 p 次，边界基函数不满足 Kronecker 性质，强施加边界条件时约束的不是真正的边界值。解决：打印U的前 p+1 个和后 p+1 个元素，确认首尾值相等且重复次数正确，再检查findSpan在 u=0 和 u=1 处的返回值是否落在合法范围。

4.2 刚度矩阵不对称或求解报奇异

现象：K组装完norm(K-K')是 1e-3 量级，或者lu分解报矩阵奇异。原因通常是单元连接关系里控制点编号重复或越界，kron顺序和dofMap不匹配，或者某个单元detJ为负。解决：先跑any(detJVec <= 0)定位坏单元，再检查elementConn每行是否有重复编号，最后核对dofMap的映射规则和B矩阵构造顺序是否一致。我一般会在组装前加一句assert(numel(unique(conn)) == numel(conn))。

4.3 边界条件施加后残差很大

现象：res在 1e-4 以上，位移场在约束点附近有明显异常。原因可能是强施加时只置了行没置列，或者置零操作把稀疏矩阵结构破坏了。解决：确认K(bcDof,:) = 0和K(:,bcDof) = 0都执行了，对角置一用speye而不是eye，右端项在置零之后赋值。如果用的是罚函数，罚因子取1e8 * max(abs(diag(K)))量级，太小约束不住，太大条件数爆炸。

4.4 后处理云图在单元边界有缝

现象：位移云图在单元交界处出现一条细线或跳变。原因是采样点恰好落在单元边界，findSpan返回的索引和相邻单元不一致，插值结果不连续。解决：采样点用linspace(0,1,nSample+2)然后去掉首尾，或者把参数点往单元内部偏移 1e-8。更稳妥的做法是按单元循环采样，每个单元内部取点，边界点用相邻单元平均值。

4.5 高次基函数下求解慢到不可接受

现象：p 从 2 加到 4，自由度没涨多少，求解时间翻十倍。原因是刚度矩阵条件数随 p 指数恶化，MATLAB 反斜杠的默认 pivoting 在病态矩阵上效率低。解决：显式用ichol做不完全 Cholesky 预条件，配合pcg，容差设 1e-10；或者对刚度矩阵做对角缩放K = D\K/D，把条件数压下来再求解。p 大于 5 时建议改用gmres带重启，别硬扛。

5. 进阶用法：把这份程序改成自己的算法试验台

跑通标准算例只是起点，这份程序真正的价值在于它结构清晰，改起来方便。我一般会做三件事：第一，把basisFuns抽成独立函数文件，换基函数（比如 T-spline 或 B 样条局部加密）时只改这一个文件，主程序不动；第二，在单元循环里加一个integrandHook函数句柄，默认是标准B'*D*B，想试新积分方案或者加材料非线性时直接替换句柄，不用改组装逻辑；第三，把postProcess里的采样和插值拆开，采样点单独存成samplePts.mat，换后处理方式时不用重算位移场。

% 可替换积分被积函数的钩子写法 integrandHook = @(B, D, detJ, w) B' * D * B * detJ * w; % 单元循环内 Ke = Ke + integrandHook(B, D, detJ, w); % 想试罚函数或者 Nitsche 时，换一个句柄即可

验证改动的正确性，我习惯用 patch test：取一个常应变场，看数值解是否精确复现。IGA 的 patch test 比有限元严格，因为基函数不是插值的，常应变场下控制点位移得按基函数值反算，不能直接给。具体做法是在边界上给线性位移场，内部控制点位移用最小二乘拟合基函数值，跑完看norm(u - u_exact)/norm(u_exact)是否在 1e-12 量级。这个测试能抓出基函数归一化错误、雅可比符号错误、组装索引错位这几类问题，比对着云图肉眼看好使。

另一个实用技巧是收敛性研究脚本化。写一个外层循环，p 从 2 到 5，每个 p 下把参数域均匀加密（h-refinement），记录自由度和能量范数误差，画 log-log 图看收敛斜率。IGA 的 k-refinement（先升次再加密）和 h-refinement 收敛率不一样，脚本化之后对比很直观。我一般会把结果存成convStudy.mat，下次改算法直接重跑对比，不用从头搭。

注意：k-refinement 时节点矢量插入要保证重复度不超过 p，升次操作后原来的节点重复度含义变了，直接复用旧矢量会出问题。常见做法是升次后重新生成 clamped 矢量，再按新次数插节点。

从那以后我每次拿到一份 IGA 代码，都先跑 patch test 再跑收敛性，两个都过了才敢往上面加新东西。这份 MATLAB 程序的好处是规模不大，改坏了容易回滚，适合当算法试验台用。希望帮到你。

本文还有配套的精品资源，点击获取