solve.h

这是等参单元的有限元程序

//////////////////////////////////////////////////////////////
// 
//    solve 是求解有限元平衡方程KV=P的函数       
//   
//    jd—是存储总体刚度矩阵主对角元在一维存储ZK数组的地址的数组
//	  zk—是一维存储的刚度矩阵，zk[0]不用来存储，赋值为零
//	Load—是存储节点等效荷载的数组，长度为NF＋1，Load[0]不用来存储，值为零
//	  NF—结构的总体自由度数
    
//////////////////////////////////////////////////////////////
#ifndef __SOLVE_
#define __SOLVE_
#include <math.h>
void Solve(int *jd,double *zk,double *Load,int NF);

void Solve(int *jd,double *zk,double *Load,int NF)

{  int    I,L,MI,LI,J,M,LJ,J2,MJ,IJ,K,J1,II,I1;
   double Z;
     for(I=1;I<=NF;I++)
	 {     L=jd[I]; 
	        if(I==1)
				MI=I+1-L;
			   else   MI=jd[I-1]+I+1-L;
           LI=L-I;
           for(J=MI;J<=I;J++)
		   {    M=jd[J];
		        LJ=LI+J;
				 J2=J-1;
                 if(J==1)
					 MJ=J+1-M;
				   else  MJ=jd[J2]+J+1-M;
                IJ=MJ; 
				if(MJ<MI) 
					IJ=MI; 
				Z=0.0;
				J1=M-J;
               if(IJ<=J2)
				   for(K=IJ;K<=J2;K++)
					   Z+=zk[K+LI]*zk[K+J1];
               Z=zk[LJ]-Z;
               if(I==J)  
				   zk[LJ]=sqrt(fabs(Z));
			        else   zk[LJ]=Z/zk[M];
		   }
           Z=0.0;
		   J1=I-1;
		   if(MI<=J1)
			   for(K=MI;K<=J1;K++)
				   Z+=zk[LI+K]*Load[K];
          Load[I]=(Load[I]-Z)/zk[L];      
	 }
     for(II=1;II<=NF;II++)
	 {   I=NF-II+1;
	     L=jd[I];
		 LI=L-I;
         if(I==1) 
			 MI=1;
		     else    MI=jd[I-1]-LI+1;
         Z=Load[I]/zk[L]; 
		 Load[I]=Z; 
		 I1=I-1;
         if(MI>I1) break;
		     else 
				 for(J=MI;J<=I1;J++)
                     Load[J]-=zk[LI+J]*Z;
	 }
}
#endif