基于哈希函数的有限元刚度矩阵组装及储存数据结构

摘要:

在有限元程序中准确快速的定位单元刚度矩阵中每一个元素在总刚矩阵中的位置对现代有限元程序运行效率的高低起着关键作用。本文主要研究有限元方法中基于稀疏结构的刚度矩阵组装问题。研究表明采用哈希函数对刚度矩阵组装比不采用该方法更能有效的组装刚度矩阵，所用的时间是不采用哈希方法所用时间的四分之一到三分之一。

关键词：

有限元刚度矩阵组装 稀疏矩阵 哈希函数

中图分类号： O174 文献标识码： A

1 引言

由于有限元方法有着广泛的工程应用，因此对有限元方法的优化在过去几十年有着深入的研究。在有限元分析中，组装单元刚度矩阵和线性方程组求解在整个分析过程中同等重要。因此对有限元方法的改进也主要集中在这两个方面，即线性方程组解法的研究和刚度矩阵的组装。线性方程组的求解涉及很多领域，其解法也日臻完善。然而刚度矩阵组装为有限元分析特有的过程，与线性方程组求解相比，其研究相对较少。因此，通过改进刚度矩阵组装方法来提高有限元程序执行效率，有很大的研究空间。本文主要研究刚度矩阵的组装，并给出基于哈希函数的刚度矩阵组装方法。

由单元刚度矩阵组装所形成的总刚矩阵是对称稀疏矩阵。[[[1] I.Z.Reguly; M.B.Giles; Finite Element Algorithms and Date Structure on Graphical Processing Units[J]. International Journal of Parallel Programming, Springer, 2013.

]]因此利用总刚矩阵的稀疏性和对称性对刚度矩阵的组装和存储方案进行优化。目前组装单元刚度矩阵主要按以下三个步骤实施。第一，选择稀疏矩阵存储方案并建立稀疏矩阵结构；第二，根据查询单元刚度矩阵中每个元素不同的指标来定位该元素在总刚矩阵的位置；第三，循环遍历每个单元，在每个单元矩阵中，把具有与总刚矩阵相同行列指标的元素加到总刚矩阵对应行列指标的位置，完成刚度矩阵的组装。尽管刚度矩阵的组装过程不算复杂，但是必须花费大量的时间去查找单元刚度矩阵中的每个元素在总刚矩阵中的对应位置，很大程度上影响有限元程序执行效率。提高查询效率可以采取提高计算机硬件性能的办法，近期有很多学者采用多核并行计算的方法来提高查询和运算效率。[[[2]Bolz,J.,Farmer,I., Grinspun,E., Schroder,P;Sparse matrix solvers on GPU:Conjugate gradients and multigrid.[J]ACM Transacations on Graphics,2003.], [[3]Komatitsch, D., Micha, D., Erlebacher, G.;Porting a high-order finite-element earthquake modeling application to NVIDA graphics cars using CUDA.[J] Journal of Parallel and Distributed Computing.2009 ], [[4] Komastisch, D., Goddeke,D.,Erlebacher,G., Michea,D.;Modeling the propagation of elastic waves using spectral elements on cluster of 192 GPUs. [J] Computer Science Research and Development,2010], [[5]Flaig, C., Arbenz,P.;A scalable memory efficient multigrid solver for micro-finite element analyses based on CT images.[J] Parallel Computing.2011]]但在这些论文中并没有涉及对刚度矩阵组装的稀疏方法的改进，它们所采用的刚度矩阵组装方法仍然是遍历行列指标查找的传统方法。再者，计算机硬件性能也不可能无限度的提高。因此，本文采用建立哈希表的方法查找单元刚度矩阵元素并确定其在总刚矩阵中的位置。通过算例验证，采用建立哈希表的方法定位单元刚度矩阵元素并组装单元刚度矩阵所花费的时间是采用循环遍历查找方法所花费时间的四分之一到三分之一。显著的提高了刚度矩阵的组装效率。

2 稀疏矩阵存储数据结构下的数据查询及单元刚度矩阵组装

在有限元分析过程中，一旦我们将结构离散化处理之后，即对单元和节点进行编号后，就确定结构的总刚矩阵结构。首先，我们选择恰当的稀疏矩阵储存方案建立稀疏矩阵储存结构，然后根据所建立的稀疏矩阵存储结构对单元刚度矩阵中的元素进行查询和计算，最后完成总体刚度矩阵的组装。常用的稀疏矩阵存储结构有坐标型方案（Coordinate,COO)、按行压缩型方案（Compressed Sparse Row, CSR)、按列压缩方案（Compressed Sparse Column, CSC) 和改进型方案（Modify Sparse Row，MSR）等。[[[6]Yousef Saad; Numberical Methods For Large Eigenvalue Problems.[M] SIAM, 2011]]为了讨论方便我们假设在结构离散化之后得到刚度矩阵K，其结构为图1。注，以下所有讨论中，编号起点均从零开始。

图1黑色方块表示需要存储数据的位置，空白处均为0。

2.1 坐标型（COO）方案

根据图1的结构我们建立三个向量AR，AC，AV来表示该总体刚度矩阵的数据结构，存储方案见表1。其中，

AR：表示矩阵K的行指标；

AC：表示矩阵K的行指标；

AV：表示矩阵K的元素。

表1

AR和AC中的值类似坐标功能，可以不用按照一定的顺序存储。采用COO型存储方案时进行单元刚度矩阵组装时，已知单元刚度矩阵某一元素在总刚矩阵的行指标和列指标需要与稀疏存储结构的行列指标进行对比，即需要遍历AR中的行指标和AC中的列指标，以确定该元素在存储结构中的位置。每一个单元刚度矩阵的每一个元素都将进行这样的操作。在实际工程中，结构离散后的单元和节点数往往数量庞大，如果再采用该数据结构进行单元刚度矩阵元素在总刚矩阵中的定位必然影响有限元程序的执行效率。因此为了提高效率，目前存储总体刚度矩阵往往采用CSR方案。

2.2 按行压缩（CRS）方案

根据图1的矩阵结构，同样我们建立三个向量AR，AR和AV。总体刚度矩阵的存储方案见表2。其中，

AR：表示矩阵K中每行第一个非零元素在这个非零元素中的编号；

AC：表示矩阵K的行指标；

AV：表示矩阵K的元素。

表2

采用CSR存储方案组装单元刚度矩阵时，在已知单元刚度矩阵中某一元素列号的时候，可以在表2的结构中快速定位出该元素在总体刚度矩阵的行指标，然后遍历该行的所有列指标与已知的列指标进行比较从而确定该元素在总刚矩阵中的位置。较COO方案中全部遍历行列指标进行比较的方法其效率有了很大提高。例如，给定单元刚度矩阵的某一元素的行指标为7，列指标为5。我们只知道7介于6和10之间，即7一定在第二行，在遍历第二行所有的列指标，我们知道5正好为第二行的第二个元素，由此便确定该元素在总体刚度矩阵的位置。 采用C++表述算法见代码1，

代码1

2.3 CSR查询方法的改进

2.3.1 哈希函数简介

哈希函数主要以哈希表的形式根据关键字对数据信息的快速查找。简而言之，哈希函数是关键字到指标的映射，指标在哈希表中作为存储数据的查找索引，通常通过数组或动态数组实现哈希表。更多关于哈希函数的介绍可以参考相关文献[[[7]Koheim, G. A. Hashing in Computer Science[M]. Canada:JOHN WILEY&SONS,INC.,

PUBLICATION,2010;]]。

2.3.2 哈希方法查询单元刚度矩阵元素

稀疏矩阵的存储方案我们仍然采用CSR方案，我们采用建立哈希函数的方法对单元刚度矩阵元素进行查找。图2为哈希函数查找单元刚度矩阵元素原理。

图2

根据图2原理，采用C++建立哈希表的算法见代码2.

代码2

在我们建立好哈希表之后，便可利用哈希表实现单元刚度矩阵元素在总刚矩阵的定位。采用C++的实现方法见代码3

代码3

2.4 单元刚度矩阵组装

在2.3中我们讨论了在稀疏矩阵存储的数据结构和元素查询，再此基础上完成刚度矩阵的组装，即把总刚矩阵中对应位置的单元刚度矩阵的元素相加，便完成单元刚度矩阵组装。采用C++实现方法见代码4。

代码4

3 数值实验

3.1 方法概述

我们采用采用六面体20节点单元对一个立方体进行离散，通过增加划分单元数目比较采用哈希表和不采用哈希表刚度矩阵的组装时间并验证采用哈希函数方法可以有效地提高刚度矩阵的组装效率。

3.2 实验计算机主要性能参数

处理器: Intel(R)Core i3-3220，3.3GHz

内存：8GB DDR3 1600MHz

操作系统： Windows7 旗舰版64位 SP1

3.3 单元划分方案及单元刚度矩阵组装

单元划分方案见表3，

表3

单元刚度矩阵组装算法采用C++实现得和核心代码，见代码5，

其中，

ne: 单元数目；

ii：总刚矩阵的行指标；

jj: 总刚矩阵的行指标；

n: 每个单元的节点数目；

m_Ks: 总刚矩阵；

ready_go: 记时起点；

finish: 记时终点；

AK.Add(ii,jj,ke[i][j]: 采用哈希表的单元刚度矩阵组装方法；

AK.Add_prime(ii,jj,ke[i][j]): 不采用哈希表的单元刚度矩阵组

装方法；

total_time: 组装单元刚度矩阵总的时间。

代码5

3.4 实验结果

通过采用在代码5中采用不同的单元刚度矩阵求和函数AK.Add(ii,jj,ke[i][j]与AK.Add_prime(ii,jj,ke[i][j])得到在不同和单元划分方案下刚度矩阵的组装计算时间见图3

图3

4 讨论与结语

通过数值实验可以发现采用建立哈希表的方法用于对单元矩阵元素的查找和刚度矩阵的组装有效地提高了有限元程序的运行效率。在图3中对实验数据的拟合处理表明应用哈希表的方法进行刚度矩阵的组装所花费的时间随着单元数的增加呈现出线性变换，其变化率约为0.，该算法的时间复杂度为O(n)。该算法以多项式时间为界限有效算法[[8]周培德 算法设计与分析[M].北京：机械工业出版社1996.2]。而采用遍历行列指标的方法时，所耗费的时间随着单元数增加的变化率为0.47。比较两种算法的变化率可以得出：采用遍历行列指标的方法所花费的时间约哈希表方法所花时间的3.3倍。

我们将哈希函数引入刚度矩阵组装为优化有限元分析程序提供了新的方法。我们相信如果将这一方法和其他方法结合，有限元程序会更加稳定高效，为我们计算更大更复杂的工程结构提供支持。

参考文献

[1] I.Z.Reguly; M.B.Giles; Finite Element Algorithms and Date Structure on Graphical Processing Units[J]. International Journal of Parallel Programming, Springer, 2013.

[2]Bolz,J.,Farmer,I., Grinspun,E., Schroder,P;Sparse matrix solvers on GPU:Conjugate gradients and multigrid.[J]ACM Transacations on Graphics,2003.

[3]Komatitsch, D., Micha, D., Erlebacher, G.;Porting a high-order finite-element earthquake modeling application to NVIDA graphics cars using CUDA.[J] Journal of Parallel and Distributed Computing.2009

[4] Komastisch, D., Goddeke,D.,Erlebacher,G., Michea,D.;Modeling the propagation of elastic waves using spectral elements on cluster of 192 GPUs. [J] Computer Science Research and Development,2010

[5]Flaig, C., Arbenz,P.;A scalable memory efficient multigrid solver for micro-finite element analyses based on CT images.[J] Parallel Computing.2011

[6]Yousef Saad; Numberical Methods For Large Eigenvalue Problems.[M] SIAM, 2011

[7]Koheim, G. A. Hashing in Computer Science[M]. Canada:JOHN WILEY&SONS,INC.,

PUBLICATION,2010;

[8]周培德 算法设计与分析[M].北京：机械工业出版社1996.2