高性能 C# 有限元

引言

开源有限元（FE）代码通常是用C++编写的，并且是带有厚厚用户手册的大型库。这并不奇怪，因为有限元计算很复杂，需要高性能。托管代码的目的不是与C在性能上竞争，因此C#之类的语言不太可能是科学计算的候选者。但有限元的情况是这样的吗？

对于热力学、电磁学、连续介质力学和其他物理学领域中的各种模型，都存在有限元公式。最常见的是，它们将线性方程组 \(\mathbf{\tilde{K} \tilde{u}=\tilde{f}}\) 组装并求解一次或多次。此步骤计算成本最高，因为矩阵 \(\mathbf{\tilde{K}}\) 通常很大。

求解线性方程组是一个算法问题，与有限元关系不大，因此该任务通常委托给单独的算法。因此，有限元代码本身可以用C#编写，利用托管集合和LINQ查询操作，最大限度地减少内存碎片。

本文描述了一种固体力学中成熟技术——*线性静态分析*的实现。当外部载荷施加到可变形物体上时，它们会在材料中产生内部应力和应变，从而产生反作用力和物体的变形。静态分析确定物体的新平衡位置和新形状。大多数商业有限元产品都具有此类功能，并解决了结构工程中的实际问题。

包含的示例网格故意做得很大，包含100,000多个单元，用于对矩阵组装过程进行基准测试。为了处理线性方程组，会调用一个第三方线性求解器，该求解器配置为使用支持CUDA的GPU或CPU。

要测试随附的演示，请按住鼠标右键并移动鼠标以变形物体。释放按钮后，求解器将立即运行。用鼠标左键旋转样本以查看结果。切换样本时，程序可能会冻结几秒钟以解析网格文件。

线性弹性有限元方法（FEM）

在有限元分析中，可变形固体由一个网格表示，该网格由节点组成，这些节点跟踪其自身相对于原始位置的位移（所谓的*拉格朗日*网格）。四面体单元允许对作用在节点上的力进行特别简单的公式化。对于每个单元，这些力等于

这里，\(\mathbf{u}=\left[u_{x1},u_{y1},u_{z1},u_{x2},u_{y2},u_{z2},u_{x3},u_{y3},u_{z3},u_{x4},u_{y4},u_{z4}\right]{}^T\) 是节点的12维位移向量，\(\mathbf{f}_{int}=\left[f_{x1},f_{y1},f_{z1},f_{x2},f_{y2},f_{z2},f_{x3},f_{y3},f_{z3},f_{x4},f_{y4},f_{z4}\right]{}^T\) 是产生的弹性力的12维向量，而 \(\mathbf{K}\) 是所谓的*单元刚度矩阵*，大小为12x12。负号强调了力相对于变形方向的方向相反（就像拉伸弹簧一样）。

对于不熟悉该主题的读者来说，刚度矩阵的大小（144个元素）可能看起来是巨大的过度。毕竟，对于弹簧模型，力和位移之间的关系要简单得多。但在固体材料中，*任何*点在*任何*方向上的**微小**变形都会在*整个*物体中产生*所有*方向的内部力。12x12的刚度矩阵捕捉了每个四面体单元的这种复杂性。这种类型的弹性响应被称为*线性弹性模型*。

这里的线性意味着力与每个单元以及整个身体中的位移呈线性比例关系。因此，每个单元的方程可以相加组合成一个描述整个对象的巨大方程组

其中 \(\mathbf{\tilde{u}}=\left[u_{x1},u_{y1},u_{z1},...,u_{zN}\right]{}^T\) 是*所有*节点（可能是数百万个）的组合位移向量，\(\mathbf{\tilde{f}}_{int}=\left[f_{x1},f_{y1},f_{z1},...,f_{zN}\right]{}^T\) 是相同大小的相应力向量，而 \(\mathbf{\tilde{K}}\) 是大小为(3N)x(3N)的*全局刚度矩阵*。一个重要的观察是，每个新节点都会向位移和力向量添加3个分量，并向全局刚度矩阵添加9个分量。

对于某些网格，矩阵 \(\mathbf{\tilde{K}}\) 非常大，以至于无法作为二维数组放入计算机内存。幸运的是，\(\mathbf{\tilde{K}}\) 是一个*稀疏*矩阵，其大多数条目为零，可以以压缩行存储（CSR）格式存储。

到目前为止，还没有提到12x12的单元刚度矩阵是如何定义或计算的。它们表示为以下乘积

其中 \(V\) 是四面体单元的体积，\(\mathbf{E}\) 是6x6的弹性矩阵（参见胡克定律），而 \(\mathbf{B}\) 是所谓的*应变-位移*矩阵，大小为6x12

系数 \(a_1,a_2,a_3,a_4,b_1,b_2,b_3,b_4,c_1,c_2,c_3,c_4\) 从四面体*雅可比*矩阵的逆中找到

其中 \(x_1,y_1,z_1,x_2,y_2,z_2,x_3,y_3,z_3,x_4,y_4,z_4\) 是四面体非位移节点的坐标。上述关系的推导可以在这些讲义中找到。

最后，以下平衡方程定义了内力和外力的平衡

其中 \(\mathbf{\tilde{f}}_{ext}\) 是外部力的向量，例如重力、摩擦力和边界上的外部载荷。在静态问题中，平衡变形是未知的。通过求解线性方程组，我们得到变形，在这种变形中，弹性响应抵消了外部载荷。力平衡表示为

关于边界条件的说明

在某些力学问题中，边界条件被表述为在物体特定部分的预设位移，例如固定在河流两侧的桥梁。在离散模型中，与固定节点对应的 \(\mathbf{\tilde{u}}\) 的分量是常数，而不是未知变量。为了处理这种情况，从线性方程组中移除固定节点，并将其对剩余自由节点的影响作为外部力来考虑。

压缩行存储（CSR）格式

如前所述，生成的矩阵太大，无法作为二维数组放入内存。CSR是存储稀疏矩阵的常用方法之一。矩阵条目由3个数组描述

public int[] rows, cols;  // structure of the sparse matrix
public double[] vals;      // non-zero values
public int N, nnz;

矩阵的大小为 N*N，*非零条目数*为 nnz。vals 数组仅包含从左到右、从上到下排列的非零元素，其大小为 nnz。cols 数组包含元素的列索引，因此其大小也等于 nnz。最后，rows 包含每行第一个非零元素的索引。rows 数组的大小为 N+1，并且 rows[N] = nnz。

在此项目中，稀疏矩阵与线性方程组的右侧（RHS）一起在 CSR_System 类中分配。rhs 向量只是一个大小为 N 的数组。所有实值都以双精度存储，但在传递给求解器之前可能转换为单精度。

矩阵组装回顾

创建全局刚度矩阵的步骤是

1. 对于每个节点，查找其邻居。
2. 仅列出非固定的、可以自由移动的节点。
3. 未知坐标的数量将是 activeNodes.Count * 3。矩阵的非零条目数将是活动节点总邻居数乘以 9。
4. 根据上一步的信息，分配CSR数组。
5. 遍历活动节点的所有邻居，在CSR矩阵的结构中创建3x3的条目。
6. 遍历所有单元，将单元刚度矩阵组装到全局刚度矩阵和RHS中。

之后，系统将进入求解器。

实现

有限元网格存储在 Model 类中的两个列表中

public List<Node> nodes;
public List<Element> elems;

节点散布在具有不同3D坐标值（x,y,z）的地方。单元是4节点四面体，节点存储在 Node[] vertices 中。在 Model() 构造函数填充 nodes 和 elems 之后，执行步骤1-5。

创建CSR结构（步骤1-5）

过程的核心是 InitializeCSR() 函数。它首先按顺序编号非固定节点

activeNodes = nodes.FindAll(nd => !nd.anchored);
int id = 0;
foreach (Node nd in activeNodes) { nd.altId = id++; }

此时，节点“不知道”它们属于哪些单元。下一段代码填充每个节点在其父单元中的列表

foreach (Element elem in elems)
    foreach (Node nd in elem.vertices) nd.AddElem(elem);

接下来，节点准备查找其邻居（通过父单元边缘连接的相邻节点）

Parallel.ForEach(nodes, nd => nd.InferConnectivityInformation());

邻居存储在 HashSet<int> 中，其中 <int> 对应于 activeNodes 中的索引。在这种情况下，整数索引比引用更好，因为它将在以后用于创建CSR结构。此步骤可以并行完成，因为没有写入冲突。

为了分配CSR，会计算活动节点总数

int count = 0;
foreach (Node nd in activeNodes) count += nd.neighbors.Count;
csr = new CSR_System(activeNodes.Count * 3, count * 9);

此时，CSR已分配，但没有结构。因此，尚无法将任何实际值写入矩阵。我们通过填充 csr.cols 和 csr.rows 来继续。此步骤包括遍历每个活动节点的每个邻居，并按顺序在CSR矩阵的结构中创建3x3条目。请注意，CreateCSRIndices(count, csr.cols) 在 csr.cols 数组中创建条目。

count = 0; // counts processed neighbors
foreach(Node nd in activeNodes) 
{
    int row_nnz = nd.CreateCSRIndices(count, csr.cols);
    csr.rows[nd.altId * 3] = count;
    csr.rows[nd.altId * 3 + 1] = count + row_nnz;
    csr.rows[nd.altId * 3 + 2] = count + row_nnz * 2;
    count += row_nnz * 3;
}

此时，CSR结构已创建，控制权将返回给GUI，用户可以在其中用鼠标修改边界条件。一旦释放鼠标右键，AssembleLinearSystem() 就会启动。

组装矩阵（步骤6）

AssembleLinearSystem() 实现上述有限元部分概述的过程。它的作用是将单元刚度矩阵分配到全局矩阵中。如果此操作并行进行，则会发生写入冲突。尽管此步骤存在并行实现，但考虑到求解系统需要更多工作，总体性能优势微不足道。

组装过程仅仅是一个嵌套循环

foreach(Element elem in elems)
{
    double[,] K = elem.K; // element stiffness matrix of size 12x12
    for(int i=0;i<4;i++) 
        for(int j=0;j<4;j++)
        {
            // disperse the entries into CSR and RHS
        }
}

一旦组装步骤完成，生成的数组将通过以下函数提交给Paralution的预条件共轭梯度求解器

        [DllImport("SolverWrapperNativeCode.dll", CallingConvention = CallingConvention.Cdecl)]
        public static extern double solveCSRDouble(int[] row_offset, int[] col, double[] val,
    int nnz, int N, double[] _rhs, double[] _x);

SolverWrapper 配置Paralution并分配适当的对象。返回值存储在 _x[] 中，然后传输到相应的节点进行渲染。

最后说明

构建项目需要以下库：Accord.NET、OpenTK和Paralution。此外，在使用CUDA设备时，还会链接CUDA Toolkit和ManagedCUDA。

此项目并非用于实际的有限元分析，但可作为编写更高级的有限元公式的起点，例如协旋模型、超弹性、隐式动力学、塑性和重网格划分。鼓励读者尝试使用使用CSR格式的其他求解器，例如PETSc。

即使对于相同的网格，求解时间也极易变化。在一台台式机上，求解器可以轻松处理具有3000万多个非零条目的矩阵，无论是CPU还是GPU。

C#中可用的多功能库在执行科学计算方面非常方便。例如，单元刚度矩阵的初始化几乎在一行代码中完成

K = B.TransposeAndMultiply(E).Multiply(B).Multiply(J.Determinant() / 6D); // K = (Bt x E x B)*V

虽然上述代码非常慢（尤其与CUDA实现相比），但对于快速无错误科学计算而言，简单性胜过效率。

源代码和演示

项目和演示可在GitHub上找到。由于其体积庞大，.geo 网格和CUDA库分别打包。只需解压所有文件即可运行演示。