HBMAT方法提取结构总刚度矩阵

phoenixfang

王新敏编的《ANSYS工程结构数值分析》中的一个例子
代码如下：

!EX6.10 提取结构刚度矩阵及处理
FINISH$/CLEAR$/FILNAME,HBFILE$/PREP7                        !定义工作文件名HBFILE.TXT
ET,1,BEAM3$MP,EX,1,2E5$R,1,1E-2,32E-5,0.5                !定义单元类型、材料、实常数
N,1$N,2,0,4$N,3,4,4$N,4,4,0$EN,1,2,3$EN,2,1,2$EN,3,4,3        !按图6-9所示创建有限元模型
F,2,FX,5$SFBEAM,1,1,PRES,10,,,,2,-1$SFBEAM,2,1,PRES,3        !施加节点荷载和单元荷载
D,1,ALL$D,4,ALL                                                !施加约束
/SOLU$SOLVE$FINISH$/AUX2                                !进入求解层求解后再进入AUX2处理器
FILE,HBFILE,FULL                                        !指定文件HBFILE.FULL
HBMAT,HBFILE,TXT,,ASCII,STIFF,YES                        !转换刚度矩阵和右边项为HBFILE.TXT文件
FINISH
!以下从HBFILE.TXT读入数据，并还原为满矩阵存储
*DIM,CONTLINE,,5                                        !定义一维数组
*VREAD,CONTLINE(1),HBFILE,TXT,,,5,,,1                        !跳过第1行后读入5个数据
(5F14.0)
PTRCRD=CONTLINE(2)                                        !保存列指针总行数
INDCRD=CONTLINE(3)                                        !保存行索引总行数
VALCRD=CONTLINE(4)                                        !保存矩阵元素总行数
RHSCRD=CONTLINE(5)                                        !保存右边项总行数
*VREAD,CONTLINE(1),HBFILE,TXT,,,4,,,2                        !跳过第2行后读入4个数据
(A3,11X,4F14.0)
NROW=CONTLINE(2)$NCOL=CONTLINE(3)                        !保存刚度矩阵的行列数
STRLINE=$CONTLINE=                                        !删除数组
*IF,RHSCRD,EQ,0,THEN                                        !如果无右边项取LS0=4行，否则取LS0=5
LS0=4$*ELSE$LS0=5$*ENDIF
*DIM,POINTR,,PTRCRD                                        !定义列指针数组
*DIM,ROWIND,,INDCRD                                        !定义行索引数组
*DIM,VALUES,,VALCRD                                        !定义矩阵元素值数组
*DIM,RHSVAL,,RHSCRD                                        !定义右边项元素值数组
*VREAD,POINTR(1),HBFILE,TXT,,,PTRCRD,,,LS0
(F14.0)                                                                 !读入列指针数据
*VREAD,ROWIND(1),HBFILE,TXT,,,INDCRD,,,LS0+PTRCRD
(F14.0)                                                                 !读入行索引数据
*VREAD,VALUES(1),HBFILE,TXT,,,VALCRD,,,LS0+PTRCRD+INDCRD
(D25.15)                                                                 !读入矩阵元素数据
*VREAD,RHSVAL(1),HBFILE,TXT,,,RHSCRD,,,LS0+PTRCRD+INDCRD+VALCRD
(D25.15)                                                                 !读入右边项元素数据
*DIM,SMATR,,NROW,NCOL                                        !定义矩阵行列数，满矩阵存储的矩阵
*DO,ICOL,1,NCOL                                                !以列数循环
STACOL=POINTR(ICOL)                                        !得到当前列指针（元素的列号）
ENDCOL=POINTR(ICOL+1)                                        !得到下一列指针
*DO,IROW,STACOL,ENDCOL-1                                !以当前列中的非零元素个数循环
TRUEROW=ROWIND(IROW)                                        !得到当前元素的行号
SMATR(TRUEROW,ICOL)=VALUES(IROW)                        !按行列号将元素值保存到矩阵中
*ENDDO$*ENDDO                                                !结束两个循环
*DO,IROW,1,NROW                                                !形成上三角元素，进而得到满矩阵
*DO,ICOL,1,NCOL
SMATR(IROW,ICOL)=SMATR(ICOL,IROW)
*ENDDO$*ENDDO
!以下为删除临时变量和数组变量
POINTR=$ROWIND=$VALUES=$RHSVAL=$ICOL=$IROW=$LS0=$STACOL=
ENDCOL=$TRUEROW=$TOTCRD=$PTRCRD=$INDCRD=$VALCRD=$RHSCRD=

用HBMAT没问题，可是运行后do loop index same as outer index ICOL，有问题，但是不会改。谁用过这个程序？

lscookie

检查下两段循环体，看是否有部分代码未识别，如*enddo，导致外层循环未结束