有限元刚度矩阵的一维变带宽存储用C++实现（一）

有限元计算中的刚度矩阵是稀疏、对称矩阵，经过各个自由度的合理排序可以使得其中的非零元素集中在主对角线附近，形成带状矩阵。这样的稀疏带状矩阵在存储时一般有“二维等带宽”和“一维变带宽”存储两种方式。其中，“一维变带宽”存储可以节省大量的内存空间，提高计算速度，但现在网上大多数“一维变带宽”存储的资料都是基于Fortran语言的，下面讨论“一维变带宽”存储的C++实现。以下仅仅是我自己学习的心得体会，如果大大们发现问题，欢迎提出讨论有限元刚度矩阵的一维变带宽存储用C++实现（一）

！

限于篇幅，关于“一维变带宽”存储的详细原理就不赘述了，有兴趣的同学可以去看有限元程序设计的书籍，下面仅讨论其C++实现。

“一维变带宽”存储需要一个存储刚度矩阵中带内元素的一维数组pGK[L](L为带内元素总数)，和一个用于存储刚度矩阵中各行主对角线上元素在数组pGK[L]中的序号（称作对角元地址）的一维辅助数组pDiag[n]（n为刚度矩阵阶数）。由于刚度矩阵中元素关于主对角线对称，故pGK[L]中只存储刚度矩阵中上三角或者下三角部分的带内元素。由于刚度矩阵中每行上非0元素仅出现在和该行所代表的自由度直接联系的那些自由度的列上，故刚度矩阵每行的最大半带宽为：

下三角存储的行最大半带宽=当前行行号-直接联系的自由度的最小列号+1；

上三角存储的行最大半带宽=直接联系的自由度的最大列号-当前行行号+1；

本篇文章仅讨论辅助数组pDiag[n]的实现，数组pGK[L]的实现将在下一篇文章中讨论，下面讲C++算法：

首先计算一维变带宽存储的辅助数组pDiag[n]，因此要想办法对节点自由度编号，之后按编好号的节点自由度组装总刚矩阵，令节点自由度编号规则为：先对未知自由度编号，再对已知自由度编号，最后把从属节点的自由度编号设置为主节点自由度编号（所谓从属节点，就是一个节点被多个单元共用的情况，就是共享节点，多个单元用的是同一个节点，自然也就用的是相同编号的自由度）。最后得到的总刚矩阵为：

有限元刚度矩阵的一维变带宽存储用C++实现（一）

也就是：下图中O代表未知位移，Θ代表已知位移，“空”代表元素为0

有限元刚度矩阵的一维变带宽存储用C++实现（一）

可以看出，指定一个当前自由度，则当前自由度的编号即为总刚矩阵中当前行的行号，当前列的列号。

得到辅助数组pDiag[n]的要点为：

（1）创建单元定位矩阵pElemDOF[nTotalElem][m](nTotalElem为单元总数；m为每个单元所拥有的自由度数，对于桁架单元有两个节点，每个节点有2个自由度，故m=4)，其中每行存储了每个单元所拥有的自由度的编号；

（2）计算下三角存储时刚度矩阵各行最大半带宽时，先取出刚度矩阵中第i行（就是第i号自由度）的编号（即，二维数组pElemDOF[i][]的第i行），之后到数组PElemDOF的第i行中寻找与该自由度相关的最小自由度编号（即，数组PElemDOF的第i行中存储的那几列上的元素的最小值），最后计算行最大半带宽。

程序流程图如下：

有限元刚度矩阵的一维变带宽存储用C++实现（一）

附代码如下：

void BandAndDiagCalcu(int nTotalElem,int nTotalDOF,Element *pElem,int ** pElemDOF,int *pDiag)

//用求出的单元定位向量矩阵求“一维压缩存储算法（一维变带宽存储）”中，刚度矩阵中各行的半带宽及对角元地址；

//最终pDiag中存储的是“对角元地址”

//求出的单元定位向量矩阵存储在pElemDOF指向的整型二维数组中
//nTotalDOF为总的自由度数；pDiag为指向整型数组的指针，该数组长度为nTotalDOF，用于存储的是“对角元地址”
//nTotalElem为总的单元数量；pElem为指向Element型的结构体数组(单元数组)的首地址
{
int iMinDOF; //节点位移(自由度)编号中的最小编号
int iBuf;
int iDOFIndex; //当前单元的当前列上的自由度编号
int i,j;
for(i=0;i<nTotalDOF;i++) //刚度矩阵中主对角线元素地址置1
pDiag[i]=1;
for(i=0;i<nTotalElem;i++) //遍历所有单元，pElemDOF[i]与(pElem+i)对应同一个单元
{
//选取当前单元的俩节点的自由度编号中的最小值放入iMinDOF
iMinDOF=pElemDOF[i][0]; //取出单元定位向量矩阵中第i行（即，第i单元）的第0号自由度
if((pElem+i)->iType==TRUSS) ////当前单元类型为桁架单元
{
for(j=0;j<4;j++) //一个桁架单元的一行在pElemDOF数组中有4列自由度
{
if(pElemDOF[i][j]<iMinDOF)
iMinDOF=pElemDOF[i][j]; //从当前桁架单元的4个节点位移编号中选出最小号
}
}
else
{
for(j=0;j<6;j++) //一个钢架单元的一行在pElemDOF数组中有6列自由度
{
if(pElemDOF[i][j]<iMinDOF)
iMinDOF=pElemDOF[i][j]; //从当前钢架单元的6个节点位移编号中选出最小号
}
}
//计算总刚矩阵中的行半带宽
if((pElem+i)->iType==TRUSS) //当前单元类型为桁架单元
{
for(j=0;j<4;j++)
{
iDOFIndex=pElemDOF[i][j]; //取出当前单元的当前列上的自由度编号
//之所以如此计算半带宽，和前面函数实现的“节点自由度编号”有关
iBuf=iDOFIndex-iMinDOF+1; //计算刚度矩阵中当前行的半带宽，
if(iBuf>pDiag[iDOFIndex]&&iBuf<=iDOFIndex+1) //判断当前节点自由度对应的半带宽。
pDiag[iDOFIndex]=iBuf;
else if(iBuf>iDOFIndex+1) //这种情况不可能发生！因为iMinDOF>=0！
pDiag[iDOFIndex]=iDOFIndex+1;
}
}
else
{
for(j=0;j<6;j++)
{
iDOFIndex=pElemDOF[i][j]; //取出当前单元的当前列上的自由度编号
//之所以如此计算半带宽，和前面函数实现的“节点自由度编号”有关
iBuf=iDOFIndex-iMinDOF+1; //计算刚度矩阵中当前行的半带宽，
if(iBuf>pDiag[iDOFIndex]&&iBuf<=iDOFIndex+1) //判断当前节点自由度对应的半带宽。
pDiag[iDOFIndex]=iBuf;
else if(iBuf>iDOFIndex+1) //这种情况不可能发生！因为iMinDOF>=0！
pDiag[iDOFIndex]=iDOFIndex+1;
}
}
}
pDiag[0]=0; //将pDiag中存储的刚度矩阵第0行的主对角元地址设为0
for(i=1;i<nTotalDOF;i++) //从pDiag数组中的第二个节点自由度编号开始循环
//pDiag[i]存储的是刚度矩阵中当前行的半带宽，pDiag[i-1]存储的是刚度矩阵中上一行的半带宽
pDiag[i]=pDiag[i]+pDiag[i-1]; //得到辅助数组
}

有限元刚度矩阵的一维变带宽存储用C++实现（一）

相关推荐