【总结】Harwell-Boeing格式矩阵转换为稀疏矩阵的方法
本帖最后由 messenger 于 2010-1-9 03:10 编辑ANSYS在8.2版及以后，增加了按照Harwell-Boeing格式输出整体矩阵（刚度、质量及阻尼矩阵）的功能。 HBMAT命令法比较简单（具体方法可参见ANSYS工程结构数值分析（王新敏，2007）一书），因而大家比较喜欢用这种方法。但是如果想将整体矩阵导入其他软件（如Matlab）进行其他运算，就有一个将Harwell-Boeing格式矩阵转换的问题。
1、Harwell-Boeing的文件格式
用HBMAT命令可输出结构刚度矩阵、质量矩阵和阻尼矩阵，其文件记录格式为大型稀疏矩阵的标准交换格式，采用索引存储方法记录矩阵的非零元素。文件基本格式是前面有4或5行描述数据，其后为单列矩阵元素值，说明如下：
第1行：格式（A72），为文件头的字符型解释，如刚度矩阵或质量矩阵等标题。
第2行：格式（5I14），分别表示该文件的总行数（不包括文件头）、矩阵列指针的总行数、矩阵行索引的总行数、矩阵元素数值的总行数、右边项总行数。
第3行：格式（A3,11X,4I14），分别为矩阵类型、矩阵行数、矩阵列数、矩阵行索引数（对组装后的矩阵，该值等于矩阵行索引数）、单元元素数（对组装后的矩阵此值为0）。
第4行：格式（2A16,2A20），分别表示列指针格式、行索引格式、系数矩阵数值格式、右边项数值格式。
第5行：格式（A3,11X,2I14），A3各列分别表示右边项格式、应用高斯起始矢量、应用eXact求解矢量；两个整数分别表示右边项列数、行索引数。三个字符中的第1个字符可取：F——全部存贮（如节点荷载向量的全部元素）、M——与系数矩阵相同方法。
第6行后：矩阵元素值（单列）。
矩阵类型用3个字符表示，第1个字符可取：R——实数矩阵、C——复数矩阵、P——仅矩阵结构（无元素数值）；第2个字符可取：S——对称矩阵、 U——不对称矩阵、H——Hermitian矩阵、Z——病态对称矩阵；R——带状矩阵；第3个字符可取：A——组装的矩阵、E——单元矩阵（未组装）。对称矩阵只存储下三角元素，如结构刚度矩阵为对称矩阵，Harwell-Boeing格式则仅记录下三角元素。
根据Harwell-Boeing文件格式，可读取矩阵的任意行列元素的数值，也可编程还原为满矩阵存储，以便它用，很显然这种提取方式比较方便。如当生成.FULL文件后，可采用命令/AUX2$FILE,mywork,full$HBMAT,mystiff,txt,ASCII,STIFF,YES$ FINISH将二进制mywork.full文件输出为ASCII码文件mystiff.txt，并输出右边项。
—— 引自 ANSYS工程结构数值分析（王新敏，2007），P357
例如：
Stiffness matrix from ANSYS FULL file dumped into Harwell-Boeing format
21 7 7 7 0
RUA 6 6 7 0
(I14) (I14) (d25.15) (d25.15)
1
3
4
5
6
7
8
1
4
2
3
4
5
6
0.106666666666667D+08
-0.533333333333333D+07
0.100000000000000D+00
0.000000000000000D+00
0.533333333333333D+07
0.200000000000000D+00
0.300000000000000D+00
其相应的含义为：
第1行为解释性文件头。
第2行的数字分别表示：该文件共有21行数据，7行列指针数据，7行矩阵行索引数据，7个矩阵元素，右边项为0行数据。
第3行：RUA表示实数矩阵、非对称、组集矩阵；该矩阵为6行6列，有7个非零元素。
第4行：指针和索引数据的输出格式均为I14，矩阵Y和右边项元素输出格式为d25.15。
第5行：因为没有右边项数据，所以没有相应的第5行。
其后为7行列指针数据、7行矩阵行索引数据，7个矩阵元素。
2、Harwell-Boeing的矩阵格式
Harwell-Boeing的矩阵的储存方法是将矩阵中的零元素去掉，然后组成一个个以列为单位的数组。
比如，上例中的矩阵的满矩阵为，
1.0667e+007 0 0 0 0 0
0 0.1 0 0 0 0
0 0 0.2 0 0 0
-5.3333e+006 0 05.3333e+006 0 0
0 0 0 0 0.3 0
0 0 0 0 0 0.4
去掉零元素后的格式为，
1.0667e+007 0.1 0.2 5.3333e+006 0.3 0.4
-5.3333e+006
要描述上面格式的数据，主要用3个量
* COLPTR —— 每一个列向量第一个元素的位置；
* ROWIND —— 每一个元素的行号；
* VALUES —— 元素值。
则，
排序位置： 1 2 3 4 5 6 7
COLPTR： 1 3 4 5 6 7 8
ROWIND： 1 4 2 3 4 5 6
VALUES：1.0667e+007 -5.3333e+006 0.1 0.2 5.3333e+006 0.3 0.4
3、Harwell-Boeing的格式的矩阵转换为Matlab的稀疏矩阵
知道了Harwell-Boeing的矩阵的储存格式以后，就可以将按照Matlab稀疏矩阵的格式转换为稀疏矩阵。
可以采用下面网址的m文件进行上述转换，http://people.sc.fsu.edu/~burkardt/m_src/hb_to_mm/hb_to_mm.m。
但这个m文件只是输出一个文本文件。
为了可以在Matlab中直接应用，对其进行了稍微修改，修改后的命令有3种形式：
kkk=hb2mm('kkk_HB.txt','kkk_mm.txt');
kkk=hb2mm('kkk_HB.txt');
hb2mm('kkk_HB.txt','kkk_mm.txt');
其中，kkk即为Matlab格式的稀疏矩阵。kkk_mm.txt为输出的含有稀疏矩阵的文本文件。kkk_HB.txt为Harwell-Boeing格式的文本文件。
再次修改：当有右边项时，原程序不能输出右边项。修改以后，当有右边项时，输出矩阵的最后一列即为右边项。
**** Hidden Message ***** 辛苦了，总结的很好 好东西啊！ 得好好研究研究！谢谢分享！！ 本帖最后由 pp786195 于 2010-1-7 16:34 编辑
再请教楼主：我用ANSYS提取的Harwell-Boeing刚度矩阵，其中输出右边项（载荷）了，但好像用MATLAB读的时候没有右边项数据？
右边项列数为1
自己试的结果：
在Harwell-Boeing文本中，MATLAB读取到数值a（刚度矩阵的最后一个数值）,则a以下的数值貌似就是我要的右边项，个数也刚好是第五行F 1**的数字。不知对否 记得当时试了试，Matlab程序可以读取到右边项。可能是你哪里用错了吧，你再检查一下吧。 fopen('8.txt') ;
K3=hb2mm('8.txt');%TXT文件在附件中
fopen('2.txt') ;
K4=hb2mm('2.txt');
用以上命令得到的K3和K4是一样的。如果有右边项应该会多一列的吧，是不是还得另外读取？我的右边项刚好就一列数值 原来原程序得到了右边项，但没有输出，现在更新了一下，当有右边项时，输出矩阵最后一列为右边项。重新下载hb2mm程序即可。 版主就是厉害，我也找到那个外国的M文件了 ，可是没发现哪段是控制右边项的。谢谢啦 本帖最后由 vibration 于 2010-2-1 19:48 编辑
你的矩阵例子解释里面可能有问题，ansys中第八个entry（COLPTR(7)=8） 应该是结束标记（其意义在于给出Scaning 的上界），你给到第七个元素（Value=0.4）上了应该是错误的。这里是HB格式的权威性说明http://people.sc.fsu.edu/~burkardt/pdf/hbsmc.pdf 有中文的，先顶在看！！！ 学习了........... 非常好，，，， 谢谢分享，这么好的东西啊。非常感谢 支持分享，谢谢 本帖最后由 darwin123 于 2011-5-6 15:59 编辑
右边项是KU=F的这个F
看了一下，将这个HARWELL-BOEING格式转化为N*N的整体刚阵还是有点复杂啊 很需要呀！谢谢了！ 正需要这个，谢谢分享 谢谢楼主，请问有人知道怎么输出特定自由度的质量矩阵的对应行吗？ 看了一下，还是不太懂啊 请问楼主COLPTR这一项的值是怎么确定的？
页:
[1]
2