求可处理3D空间中2D流形三角网格的Laplace-Beltrami特征值问题的有限元求解器

问题背景 #

我尝试使用有限元法求解嵌入在3维空间中的2维流形（例如球面或环面的边界）上的Laplace-Beltrami特征值问题。已通过Pyvista库生成三角网格，数据如下：

网格点数据 #

points = [[-2.5819006 -6.5754533 -2.0017405]
 [-2.6988158 -6.5754533 -1.9960535]
 [-2.5819006 -6.604572  -1.9960535]
 ...
 [ 1.4294114  7.86527    1.7230791]
 [ 2.2316737  7.86527    1.6120211]
 [ 3.0339363  7.86527    1.388944 ]]

网格单元数据 #

cells = [[   0    1    2]
 [   2    3    0]
 [   4    5    6]
 ...
 [ 925 1144  924]
 [1145 1144  925]
 [ 909 1145  925]]

日常使用scikit-fem处理平面2维三角网格或3维四面体网格效果良好，但该库无法直接处理3维空间中的2维曲面三角网格，现寻求可轻松解决此类网格Laplace-Beltrami特征值问题的工具包。

推荐工具包及实现思路 #

1. FEniCSx（FEniCS新一代版本） #

FEniCSx原生支持嵌入高维空间的低维流形网格，是处理这类问题的首选工具之一。

核心步骤：

• 将PyVista生成的网格转换为FEniCSx可识别的格式：

  from dolfinx import mesh
  import numpy as np
  from mpi4py import MPI
  
  # 转换单元格式：每个三角形单元前添加3（表示三角形单元的节点数）
  cells_fenics = np.hstack([np.full((cells.shape[0], 1), 3), cells])
  # 创建FEniCSx曲面网格
  domain = mesh.create_mesh(MPI.COMM_WORLD, cells_fenics, points, mesh.CellType.triangle)
  

• 定义连续伽辽金函数空间（如CG1）：

  from dolfinx.fem import FunctionSpace
  V = FunctionSpace(domain, ("CG", 1))
  

• 组装Laplace-Beltrami算子的刚度矩阵与质量矩阵：

  from dolfinx.fem import form
  from ufl import TrialFunction, TestFunction, inner, grad, dx
  
  u = TrialFunction(V)
  v = TestFunction(V)
  # 刚度矩阵对应Laplace-Beltrami算子
  a = form(inner(grad(u), grad(v)) * dx)
  # 质量矩阵用于特征值问题
  m = form(inner(u, v) * dx)
  

• 转换为稀疏矩阵并求解特征值问题（可配合SciPy或FEniCSx自带的求解器）：

  from dolfinx.fem.petsc import assemble_matrix
  from scipy.sparse.linalg import eigs
  
  A = assemble_matrix(a).toarray()
  M = assemble_matrix(m).toarray()
  # 求解前k个最小特征值
  eigenvalues, eigenvectors = eigs(A, k=10, M=M, which='SM')
  

2. PyMeshLab + SciPy #

若偏好轻量工具链，可借助PyMeshLab处理网格拓扑，手动组装Laplace-Beltrami的离散矩阵：

• 用PyMeshLab导入PyVista网格，计算每个三角形的第一基本形式（度量张量）；
• 基于度量张量计算单元刚度矩阵和质量矩阵，全局组装后用SciPy的特征值求解器计算结果。

3. Trilinos #

针对大规模高性能计算场景，Trilinos框架的Belos特征值求解器可高效处理此类问题，配合MeshKit完成曲面网格的预处理，但学习曲线相对陡峭。

补充：scikit-fem的变通方案 #

若坚持使用scikit-fem，可手动计算曲面单元的雅可比行列式（基于第一基本形式），修改单元积分逻辑。但该方案需要自定义单元类型，不如直接使用原生支持曲面的框架便捷。

内容的提问来源于stack exchange，提问作者SebastianP