testfem.cpp

计算力学编程,有限元

#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include "function.h"

void main(void)
{ char *InFile;                                             //有限元输入数据文件
  char *OutFile;                                            //有限元输出数据文件
  int NodeNum;                                                        //节点总数
  int EleNum;                                                         //单元总数
  int Iopt;                         //Iopt=0，平面应力；=1，平面应变；=2，轴对称
  int ND;                                                 //单个节点的自由度个数
  int nGauss;                                                   //高斯积分点个数
  double Thick;                                       //单元厚度，仅用于平面应力
  double E,u;                                             //单元弹性模量和泊松比
  double Density;                                                     //单元密度
  double *X,*Y;                                                   //节点坐标数组
  int *JD;                                //总刚度矩阵主对角元在一维存储中的位置
  int (*Ele)[8];                                              //单元节点编号数组
  double *zk;                                             //一维存储的总刚度矩阵
  double *F;                                                      //节点载荷向量
  int NLoad;                                                      //集中载荷个数
  int SLoad;                                                      //压力边界个数
  int (*SLoadSide)[3];                                      //压力边界的节点编号
  double (*SLoadPo)[3];                                     //压力边界的压力数据
  int MLoad;                             //=0，体积不考虑体积力；!=0，考虑体积力
  double ax,ay;                                                     //惯性加速度
  int FixedNum;                                                       //约束个数
  int (*FixedDOF)[2];                                       //被约束的自由度信息
  double *FixedData;                                                    //约束值
  double (*stress)[5];

  int i,j,k,Je[16];                                                   //临时变量
  double xo[8],yo[8],ek[256],uo[8],vo[8];                             //临时变量

  InFile="InFile.dat";
  ifstream In(InFile);                    //打开输入文件，输入文件名为InFile.dat
  if(!In)
  { ofstream Error("file.err");
    Error<<"no input file ";Error.close();exit(1);
  }                                                   //未发现输入文件infile.dat
  OutFile="OutFile.dat";
  ofstream Out(OutFile);Out.setf(ios::scientific);
  Out.setf(ios::right); Out.precision(7);Out.fill(' ');

  In>>NodeNum>>ND;                        //输入节点总数以及单个节点的自由度个数
  In>>Iopt>>nGauss;
  if(!(Iopt==0||Iopt==1||Iopt==2)) Iopt=0;                //设定单元选项的缺省值
  if(nGauss!=2) nGauss=3;                                 //设定高斯积分的缺省值
  X=new double [NodeNum+1];               //为节点X坐标分配内存空间，下标从1开始
  Y=new double [NodeNum+1];               //为节点Y坐标分配内存空间，下标从1开始
  JD=new int [NodeNum*ND+1];                     //为JD分配内存空间，下标从1开始
  F=new double [NodeNum*ND+1];                    //为F分配内存空间，下标从1开始
  stress=new double[NodeNum+1][5];           //为stress分配内存空间，下标从1开始

  for(i=1;i<=NodeNum;i++) {In>>j;In>>X[j]>>Y[j];}   //输入所有节点的坐标
  

  Out<<"\n                有限元分析数据         \n";
  Out<<"\n                  节点坐标\n";
  Out<<"  节点号     x坐标           y坐标\n";
  for(i=1;i<=NodeNum;i++)
  Out<<setw(6)<<i<<setw(16)<<X[i]<<setw(16)<<Y[i]<<"\n";

  for(i=0;i<=NodeNum*ND;i++) JD[i]=0;                               //为JD赋初值
  for(i=0;i<=NodeNum*ND;i++) F[i]=0;                                 //为F赋初值
  for(i=0;i<=NodeNum;i++) for(j=0;j<5;j++) stress[i][j]=0.0;    //为stress赋初值

  In>>EleNum;                                                     //输入单元个数
  Ele=new int [EleNum][8];                          //为单元节点编号分配内存空间

  for(i=0;i<EleNum;i++) for(j=0;j<8;j++) In>>Ele[i][j];       //输入单元节点编号
  Out<<"\n有限元单元数据\n";
  Out<<"    节点1  节点2  节点3  节点4  节点5  节点6  节点7  节点8\n";
  for(i=0;i<EleNum;i++)
  { for(j=0;j<8;j++) Out<<setw(7)<<Ele[i][j];Out<<"\n";}

  In>>Thick>>E>>u;                                  //输入厚度、弹性模量和泊松比
  if(Iopt) Thick=1.0;       

  if(Iopt==0) Out<<"\n单元类型=平面应力\n";
  if(Iopt==1) Out<<"\n单元类型=平面应变\n";
  if(Iopt==2) Out<<"\n单元类型=轴对称\n";
  Out<<"\n单元材料性质参数\n";
  Out<<"     杨氏模量      泊松比       厚度(用于平面应力)\n" ;
  Out<<setw(14)<<E<<setw(14)<<u<<setw(14)<<Thick<<"\n";
  Out<<"\n高斯积分点个数="<<nGauss<<"\n";

  In>>NLoad;                                                  //输入集中载荷个数
  if(NLoad>0)
  { Out<<"\n集中载荷数据\n";
    Out<<"节点号 自由度号      载荷值\n";
    for(i=0;i<NLoad;i++)
    { In>>j>>k; In>>F[(j-1)*ND+k];
      Out<<setw(6)<<j<<setw(6)<<k<<setw(18)<<F[(j-1)*ND+k]<<"\n";
    }
  }
  In>>SLoad;                                                 //输入压力边界的个数
  if(SLoad>0)
  { SLoadSide=new int [SLoad][3];
    SLoadPo=new double [SLoad][3];
    for(i=0;i<SLoad;i++)
    { In>>SLoadSide[i][0]>>SLoadSide[i][1]>>SLoadSide[i][2];
      In>>SLoadPo[i][0]>>SLoadPo[i][1]>>SLoadPo[i][2];
    }
    Out<<"\n表面载荷数据\n";
    Out<<"节点号 节点号 节点号      载荷1         载荷2           载荷3\n";
    for(i=0;i<SLoad;i++)
    { for(j=0;j<3;j++)  Out<<setw(6)<<SLoadSide[i][j];
      for(j=0;j<3;j++)  Out<<setw(16)<<SLoadPo[i][j];
      Out<<"\n";
    }
  }

  In>>MLoad;                                                    //输入体积力标志
  if(MLoad)
  { In>>Density>>ax>>ay;
    Out<<"\n体积力载荷数据\n";
    Out<<"      密度          x向加速度      y向加速度\n";
    Out<<setw(14)<<Density<<setw(16)<<ax<<setw(16)<<ay<<"\n";
  }

  In>>FixedNum;                                               //输入位移约束个数
  if(FixedNum<=0)
  { ofstream Error("file.err");
    Error<<"not fixed\n";Error.close();exit(1);
  }                                                   
  FixedDOF=new int [FixedNum][2];
  FixedData=new double [FixedNum];

  for(i=0;i<FixedNum;i++)
  { In>>FixedDOF[i][0]>>FixedDOF[i][1]; In>>FixedData[i];}
  Out<<"\n指定位移数据\n";
  Out<<"节点号 自由度号      约束值\n";
  for(i=0;i<FixedNum;i++)
  Out<<setw(6)<<FixedDOF[i][0]<<setw(6)<<FixedDOF[i][1]<<setw(18)<<FixedData[i]<<"\n";
  In.close();

//完成所有数据的输入
  skdd(JD,Ele,NodeNum,ND,EleNum);       //计算刚度矩阵主对角元在一维存储中的地址
  zk=new double [JD[NodeNum*ND]+1];                     //为刚度矩阵分配存储空间
  for(i=0;i<=JD[NodeNum*ND];i++)
  zk[i]=0.0;                                                    //刚度矩阵初始化

  for(i=0;i<EleNum;i++)
  { for(j=0;j<8;j++) if(Ele[i][j])
    {xo[j]=X[Ele[i][j]];yo[j]=Y[Ele[i][j]];
     Je[2*j]=2*Ele[i][j]-1;Je[2*j+1]=2*Ele[i][j];
    }                                         //返回当前单元的坐标以及自由度编号
    else {xo[j]=0;yo[j]=0;Je[2*j]=0;Je[2*j+1]=0;}
    if(Ele[i][7]>=0)
    { CacuEk(E,u,Thick,xo,yo,ek,Ele[i],Iopt,nGauss);    //计算四边形单元刚度矩阵
      Assemble(ek,Je,16,zk,JD);                               //组装单元刚度矩阵
      if(MLoad)
      CacuBF(Density,Thick,ax,ay,xo,yo,F,Ele[i],Iopt,nGauss);   //计算体积力载荷
    }
    else
    { CacuEk(E,u,Thick,xo,yo,ek,Ele[i],Iopt);           //计算三角形单元刚度矩阵
      Assemble(ek,Je,12,zk,JD);
      if(MLoad)
      CacuBF(Density,Thick,ax,ay,xo,yo,F,Ele[i],Iopt);
    }
  }                               //完成单元刚度矩阵的计算、组装以及体积力的计算

  for(i=0;i<SLoad;i++)
  { for(j=0;j<3;j++) if(SLoadSide[i][j])
    {xo[j]=X[SLoadSide[i][j]];yo[j]=Y[SLoadSide[i][j]];}  //返回当前单元边界坐标
    CacuSF(Thick,SLoadPo[i],xo,yo,F,SLoadSide[i],Iopt,nGauss);  //计算表面力载荷
  }                                                     //完成表面分布载荷的处理
  for(i=0;i<FixedNum;i++)                                         //处理边界条件
  Fix((FixedDOF[i][0]-1)*ND+FixedDOF[i][1],FixedData[i],zk,F,JD,NodeNum*ND);

  LLT(JD,zk,F,NodeNum*ND);                              //平方根法求解线性方程组

  Out<<"\n                有限元分析数据         \n";
  Out<<"\n                  节点位移\n";
  Out<<"  节点号     x位移           y位移\n";
  for(i=1;i<=NodeNum;i++)
  Out<<setw(6)<<i<<setw(16)<<F[2*i-1]<<setw(16)<<F[2*i]<<"\n";

  for(i=0;i<EleNum;i++)
  { for(j=0;j<8;j++) if(Ele[i][j])
    {xo[j]=X[Ele[i][j]];yo[j]=Y[Ele[i][j]];
     Je[2*j]=2*Ele[i][j]-1; Je[2*j+1]=2*Ele[i][j];
     uo[j]=F[Je[2*j]] ;vo[j]=F[Je[2*j+1]];
    }                                          //返回当前单元的坐标以及自由度编号
    else
    { xo[j]=0;yo[j]=0;Je[2*j]=0; Je[2*j+1]=0;uo[j]=0;vo[j]=0; }
    if(Ele[i][7]>=0)
    CacuSS(E,u,xo,yo,uo,vo,stress,Ele[i],Iopt);   //按单元计算四边形单元节点应力
    else
    CacuSS(E,u,xo,yo,uo,vo,stress,Ele[i],&Iopt);  //按单元计算三角形单元节点应力
  }
  for(i=1;i<NodeNum;i++) for(j=1;j<5;j++) stress[i][j]/=stress[i][0];  //应力平均                                                 //改过

  Out<<"\n                   节点应力              \n";
  if(Iopt)
  Out<<"\n 节点号    x向应力     y向应力     xy向应力    z向应力    \n";
  else
  Out<<"\n 节点号    x向应力     y向应力     xy向应力    \n";
  for(i=1;i<=NodeNum;i++)
  { Out<<setw(6)<<i;
    if(Iopt) for(j=1;j<5;j++) Out<<setw(16)<<stress[i][j];
    else   for(j=1;j<4;j++)   Out<<setw(16)<<stress[i][j];
    Out<<"\n";
  }

  Out.close();
  InFile=NULL;OutFile=NULL;
  if(X) delete [] X;
  if(Y) delete [] Y;
  if(JD) delete [] JD;
  if(Ele) delete [] Ele;
  if(zk) delete [] zk;
  if(F)delete [] F;
  if(stress)delete [] stress;
  if(SLoad>0) { delete [] SLoadSide; delete [] SLoadPo;}
  if(FixedDOF)  {delete [] FixedDOF;delete [] FixedData;}
  return ;
}
//------------------------