skdd.h

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

///////////////////////////////////////////////////////////////////////
//
//         skdd函数：根据节点和单元信息计算总体刚度矩阵中非零元素的排布信息，
//         Getstiff,Setstiff函数：根据二维刚度矩阵的行，列号来取对应的一维刚度矩阵中的值
//         Skke函数：  单元刚度矩阵组装程序
//         jd[i]     --表示二维总刚中第i行第i列的元素在一维总刚数组中的序号，
//                        从而建立二维总刚和一维总刚的联系。
//         Ele       --EleNumber*8的存储单元的节点编号的二维数组指针,Ele[i][j]第i行单元的第j个节点编号.
//         NodeNumber--节点总数;
//         ND        --每个节点的自由度个数,对所有节点均相同;
//         EleNumber --单元总数;
//         zk        --因为存储的刚度矩阵长度为jd[jd[0]]+1     zk[0]未用
//         ek        --存放一维按行存储单元刚度数组,长度为16*16;
//         je        --存放单元自由度总体编号数组,长度为 16
//                      如果je[i]=0,这表明单元中不存在第i个自由度
//         NFE       --单元总自由度数
//
////////////////////////////////////////////////////////////////////////////////////////
#ifndef __SKDD_
#define __SKDD_
void skdd(int *jd,int (*Ele)[8],int NodeNumber,int ND,int EleNumber);
void SetStiff(int i,int j,double Value,double *zk,int *jd);
double GetStiff(int i, int j,double *zk,int *jd);
void Skke(double *ek,int *je,int NFE,double *zk,int *jd);

void skdd(int *jd,int (*Ele)[8],int NodeNumber,int ND,int EleNumber)

{  int i,j,k,m,d,jew[8],*Min,*Band;    //定义临时变量
   Min=new int[NodeNumber+1];
   for(i=1;i<=NodeNumber;i++) 
	   Min[i]=NodeNumber;          //Min存放每个节点相连的最小节点号,此处赋初值;
   for(i=0;i<EleNumber;i++)
  {   for(j=0;j<8;j++)  
           jew[j]=Ele[i][j];
       k=jew[0];  
	   for(j=1;j<8;j++)
         if(jew[j]>0) 
			 if(jew[j]<k)
				  k=jew[j];  //寻找第i个单元中的最小节点号;
      for(j=0;j<8;j++)
      {   m=jew[j];
	      if(m>0) 
			  if(k<Min[m])  Min[m]=k;
	  }                            //更新第i个单元内所有节点的相连的最小节点号;
  }                                     //完成对单元的循环,寻找最小节点号;
  Band=new int[ND*NodeNumber+1];                  //为半带宽分配存储空间;

  for(i=1;i<=NodeNumber;i++)                      //完成所有半带宽的计算;
   {
      d=i-Min[i];
     for(j=1;j<=ND;j++)  
		 Band[ND*(i-1)+j]=ND*d+j;
   }
  jd[1]=1;  
  for(i=2;i<=ND*NodeNumber;i++) 
	  jd[i]=jd[i-1]+Band[i];
 jd[0]=0 ;  
 delete[]Min;  
 delete[]Band;
 return;
}

void SetStiff(int i,int j,double Value,double *zk,int *jd)
/* i-要访问的刚度矩阵元素的行号
   j-要访问的刚度矩阵元素的列号
   Value-刚度矩阵元素K[i][j]
*/
{     if(i<j) return;                       //上三角部分 ,不作处理
        if(j<i-jd[i]+jd[i-1]+1) return;        //带宽以外,不作处理
      zk[jd[i]+j-i]=Value;
}  
                                   //设定矩阵元素k(i,j)=Value
double GetStiff(int i, int j,double *zk,int *jd)
{   int ii,jj;
    if(i<j) 
	{  ii=j;  jj=i;} 
	   else { ii=i; jj=j;}    //上三角部分,交换下标
    if(jj<ii-jd[ii]+jd[ii-1]+1) 
		return (0);          //在半带宽以外,返回零值
       else 
		   return (zk[jd[ii]+jj-ii]);
}                                        //返回矩阵元素k(i,j)


//单元刚度矩阵组装程序

void Skke(double *ek,int *je,int NFE,double *zk,int *jd)
 {    
	 int    i,j;
       double Value;
       for(i=0;i<NFE;i++) 
		   if (je[i])  
			   for(j=0;j<NFE;j++)
                  if(je[j])
				  {   Value=GetStiff(je[i],je[j],zk,jd)+ek[i*16+j];
                      SetStiff(je[i],je[j],Value,zk,jd);
				  }  
 }



#endif