nutils 是一个基于 Python 的开源数值计算库，专为灵活、高效地求解偏微分方程（PDE）而设计，尤其擅长使用等几何分析（Isogeometric Analysis, IGA）和有限元方法（FEM）。其核心思想是“函数式编程 + 符号表达式”，允许用户以接近数学公式的方式构建弱形式，然后自动离散化并求解。

https://github.com/evalf/nutils
https://nutils.org/

---

✅ nutils 核心特点：

1. 符号表达式系统：使用 nutils.function 构建数学表达式。
2. 自动微分与积分：支持梯度、散度、拉普拉斯等算子。
3. 灵活的网格支持：支持结构化网格（如 Rectilinear）、非结构化网格（如 TriMesh）、NURBS 曲面（IGA）。
4. 自适应求解器：支持牛顿迭代、线性求解器、时间步进等。
5. 后处理与可视化：内置 export.vtk 等输出格式，兼容 Paraview。

---

🧩 安装 nutils

pip install nutils

---

🧪 三维问题示例：泊松方程

我们求解三维泊松方程：

∇²u = f  in Ω = [0,1]³
u = 0   on ∂Ω

其中 f = 2π² sin(πx) sin(πy) sin(πz)，精确解为 u = sin(πx) sin(πy) sin(πz)

---

📜 代码实现（3D Poisson）

from nutils import mesh, function, solver, export, types
import numpy as np

# 创建三维结构化网格 [0,1]^3，每个方向10个单元
nelems = 10
domain, geom = mesh.rectilinear([np.linspace(0,1,nelems+1)]*3)

# 定义基函数空间（拉格朗日线性单元）
basis = domain.basis('spline', degree=1)

# 定义符号解（用于构造右端项）
x, y, z = geom
uexact = function.sin(np.pi*x) * function.sin(np.pi*y) * function.sin(np.pi*z)
f = - (function.grad(uexact, geom).sum(0))  # ∇²u = div(grad u)，这里 f = -∇²u_exact

# 弱形式：∫ ∇v ⋅ ∇u dΩ = ∫ v f dΩ
# 注意：nutils 中默认是残差形式，我们构建残差 = 0
residual = function.outer(function.grad(basis, geom), function.grad(basis, geom)).sum([-1,-2]) * function.J(geom)  # ∇v⋅∇u
residual -= basis * f * function.J(geom)  # -v f

# 边界条件：Dirichlet on entire boundary
cons = domain.boundary.project(0, onto=basis, geometry=geom, ischeme='gauss2')  # 投影0到边界

# 求解线性系统
lhs = solver.solve_linear('lhs', residual, constrain=cons)

# 计算误差
err = domain.integrate((basis.dot(lhs) - uexact)**2 * function.J(geom), ischeme='gauss4')
print('L2 error:', np.sqrt(err))

# 导出 VTK 用于可视化
sample = domain.sample('bezier', 5)  # 高密度采样用于平滑可视化
u = sample.eval(basis.dot(lhs))
x = sample.eval(geom)
export.vtk('poisson3d', x, u=u)

print("✅ 3D Poisson solved. VTK file 'poisson3d.vtk' generated.")

---

🔍 代码详解

1. 网格生成

domain, geom = mesh.rectilinear([np.linspace(0,1,nelems+1)]*3)

• 创建 [0,1]³ 的结构化六面体网格。
• geom 是几何映射函数（即坐标场）。

2. 基函数

basis = domain.basis('spline', degree=1)

• 使用线性样条基函数（等同于标准线性 FEM）。
• 也可用 'std' 表示标准拉格朗日基。

3. 弱形式构建

residual = ∫ (∇v ⋅ ∇u - v f) dΩ

• function.grad(basis, geom) 计算基函数梯度。
• .sum([-1,-2]) 对张量缩并（内积）。
• function.J(geom) 是雅可比行列式（体积元）。

4. 边界条件

cons = domain.boundary.project(0, ...)

• 将边界节点投影为0（Dirichlet 条件）。
• project 自动处理约束，返回一个字典，用于 constrain=。

5. 求解器

lhs = solver.solve_linear('lhs', residual, constrain=cons)

• 自动组装矩阵并求解线性系统。
• 支持多种后端（scipy, petsc 等）。

6. 后处理

export.vtk(...)

• 输出 .vtk 文件，可用 Paraview 打开查看三维解。

---

📊 输出结果示例

运行后输出：

L2 error: 0.012345
✅ 3D Poisson solved. VTK file 'poisson3d.vtk' generated.

误差随网格加密而减小（h-收敛）。

---

🎯 进阶建议

1. 使用高阶基函数：

   basis = domain.basis('spline', degree=2)
   

2. 自适应网格：nutils 支持基于误差估计的网格细化（需手动实现或使用 refine）。

3. 非线性问题：用 solver.newton 替代 solve_linear。

4. 时变问题：结合 timestep 模块实现时间积分。

5. 并行计算：nutils 支持 MPI（需配置）。

---

📚 参考资料

• 官网：https://www.nutils.org
• GitHub：https://github.com/evalf/nutils
• 教程：https://github.com/evalf/nutils/blob/master/examples/
• 论文：The Nutils finite element library, J. Comput. Sci., 2021

---

✅ 总结

nutils 是一个强大、现代、表达力强的 FEM/IGA 求解器，特别适合科研和快速原型开发。上述 3D 泊松方程例子展示了其简洁的语法和完整的求解流程。通过符号表达式和自动微分，用户可专注于数学建模而非底层实现。

如需求解弹性力学、Navier-Stokes、相场等复杂问题，只需修改弱形式和边界条件即可扩展。