ADINA壳单元刚度矩阵导出

从ADINA中导出壳单元刚度矩阵，结果按照ADINA特定规范排列，通过MATLAB编写代码从ADINA导出的刚度矩阵文件读取并还原刚度矩阵

***GLOBAL STIFFNESS MATRIX, DIRECT MATRIX INPUT FILE***
dmig-define  stif-C:\Users\chenh\Desktop\ADINA\two_test_1
1	X-T	0.0000000000000E+00	col
1	X-T	0.1176271906564E+10
1	Y-T	0.0000000000000E+00	col
1	Y-T	0.8191472696014E+09
1	X-T	-0.1785623184815E+09
1	Z-T	0.0000000000000E+00	col
1	Z-T	0.8191472696014E+09
1	Y-T	-0.1190415456543E+09
1	X-T	0.1785623184815E+09
......

使用如下代码还原刚度矩阵

argv="plane.mtx";               % 文件名
infile=fopen(argv);
k=zeros(441*6);                 % 手动定义刚度矩阵大小
while ~feof(infile)
    tline=fgetl(infile);        % 按行读文件
    a=textscan(tline,'%s');
    if size(a{1},1)==4          % 确定当前刚度矩阵元素行数i
        dof=a{1}(2);
        dof=char(dof);
        switch dof
            case 'X-T'
                dof=1;
            case 'Y-T'
                dof=2;
            case 'Z-T'
                dof=3;
            case 'X-R'
                dof=4;
            case 'Y-R'
                dof=5;
            case 'Z-R'
                dof=6;
        end
        indice=a{1}(1);indice=str2double(indice);
        i=indice*6-6+dof;
    elseif size(a{1},1)==3      % 确定当前刚度矩阵元素列数j
        dof=a{1}(2);
        dof=char(dof);
        switch dof
            case 'X-T'
                dof=1;
            case 'Y-T'
                dof=2;
            case 'Z-T'
                dof=3;
            case 'X-R'
                dof=4;
            case 'Y-R'
                dof=5;
            case 'Z-R'
                dof=6;
        end
        indice=a{1}(1);indice=str2double(indice);
        j=indice*6-6+dof;
        k(i,j)=str2double(a{1}(3));     % 将当前刚度矩阵元素储存进刚度矩阵
    end
end
k=tril(k,-1)+tril(k)';          % 下三角矩阵还原为整个刚度矩阵

运行程序，得到还原刚度矩阵如下：

1.176E9	-1.179E8	1.179E8	-3.670E-11	-4.474E-10	...
-1.179E8	8.192E8	-1.190E8	5.834E7	8.251E7	...
1.179E8	-1.190E8	8.192E8	5.834E7	8.251E7	...
-3.670E-11	5.834E7	5.834E7	2.918E7	2.104E3	...
-4.474E-10	8.251E7	8.251E7	2.104E3	2.918E7	...
...	...	...	...	...	...