if(N1==0) nd=3*np; else nd=6*np; kd=new long[nd+1]; for(j=0;j<=nd;j++) kd[j]=0;
//确定每一自由度的最大半带宽及一维存储数组的最大值
for(j=1;j<=np;j++) { id1=0; for(i=0;i
if(j1!=j&&j2!=j) goto pas200; if(j1>j||j2>j) goto pas200; id2=abs(j2-j1); if(id2>id1) id1=id2; pas200:; }
//杆单元
if(N1==0) { id=3*(id1+1); if(j==1) kd[1]=1; if(j!=1) kd[3*j-2]=kd[3*j-3]+id-2; kd[3*j-1]=kd[3*j-2]+id-1; kd[3*j]=kd[3*j-1]+id; }
//定义一维存储矩阵数组
k=new float[kd[nd]+1]; kp=new float[kd[nd]+1]; for(i=1;i<=kd[nd];i++) { k[i]=0; }
//杆单元,求单刚，并组装总刚
if(N1==0) { for(ie=0;ie
inn=nln[ie]-(nln[ie]/50)*50-1; j1=nl1[ie]-1; j2=nl2[ie]-1; x1=npx[j1]; y1=npy[j1]; z1=npz[j1]; x2=npx[j2]; y2=npy[j2]; z2=npz[j2];
xx=x2-x1; yy=y2-y1; zz=z2-z1; ole=sqrt(xx*xx+yy*yy+zz*zz); s=sqrt(xx*xx+yy*yy);
ke[1][1]=nje*njA[nln[ie]-(nln[ie]/50)*50-1]/ole; ke[4][4]=ke[1][1]; ke[4][1]=-ke[1][1]; ke[1][4]=-ke[1][1]; if(s!=0) { t[1][1]=xx/ole; t[1][2]=yy/ole; t[1][3]=zz/ole; t[2][1]=-yy/s; t[2][2]=xx/s;
t[2][3]=0; t[3][1]=-t[2][2]*t[1][3]; t[3][2]=t[2][1]*t[1][3]; t[3][3]=s/ole; }
else { t[1][1]=0; t[1][2]=0; t[1][3]=1; t[2][1]=0; t[2][2]=1; t[2][3]=0; t[3][1]=-1; t[3][2]=0; t[3][3]=0; }
for(i=4;i
for(i=1;i
for(i=1;i
for(i=1;i
for(j=1;j
if(nn>ll) goto pas260; ij=kd[ll]+nn-ll; k[ij]=k[ij]+ke[n][l];
pas260:; }} }} }
这个程序得到一维变带宽存储的刚度矩阵K。如何实现K与位移向量相乘并得出计算结果
