main.cpp

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

////////////////////////////////////////////////////////////////////////////// 
//
//  程序简介：该程序是针对二维平面或轴对称问题编制的有限元程序，采用4－8节点
//
//            等参数单元，只能处理各向同性材料的简单弹性问题。
//
//  以下是主体函数，包括“原始数据文件的读入”，“主体计算”                
//////////////////////////////////////////////////////////////////////////////
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <stdlib.h>
#include <math.h>
#include "skdd.h"
#include "cacuek.h"
#include "cacubf.h"
#include "cacusf.h"
#include "fixed.h"
#include "solve.h"
#include "cacuss.h"

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

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

	       ///////////////////////////////
           // 原始数据文件的读入，分配内存
	       ///////////////////////////////
  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>>NodeNumber>>ND;                      //输入节点总数以及单个节点的自由度个数
  In>>IK>>nGauss;
  if(!(IK==0||IK==1||IK==2) )  IK=0;                   //设定单元选项的缺省值
  if(nGauss!=2) nGauss=3;                                 //设定高斯积分的缺省值

  X=new double [NodeNumber+1];               //为节点X坐标分配内存空间，下标从1开始
  Y=new double [NodeNumber+1];               //为节点Y坐标分配内存空间，下标从1开始
  JD=new int [NodeNumber*ND+1];                     //为JD分配内存空间，下标从1开始
  F=new double [NodeNumber*ND+1];                    //为F分配内存空间，下标从1开始
  Stress=new double[NodeNumber+1][5];             //为stress分配内存空间

  for(i=1;i<=NodeNumber;i++)
  {In>>j;In>>X[j]>>Y[j];}                //输入所有节点的坐标，不一定按顺序输入

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

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


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

  for(i=0;i<EleNumber;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<EleNumber;i++)
  { for(j=0;j<8;j++) Out<<setw(7)<<Ele[i][j];Out<<"\n";}

  In>>tt>>E>>u;
  if(IK) tt=1.0;               //在非平面应力情况下，将单元厚度设为单位厚度

  if(IK==0) Out<<"\n单元类型是平面应力\n";
  if(IK==1) Out<<"\n单元类型是平面应变\n";
  if(IK==2) Out<<"\n单元类型是轴对称\n";
  Out<<"\n单元材料性质参数\n";
  Out<<"     杨氏弹性模量      泊松比          厚度(仅用于平面应力)\n" ;
  Out<<setw(17)<<E<<setw(17)<<u<<setw(17)<<tt<<"\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(3)<<j<<setw(3)<<k<<setw(16)<<F[(j-1)*ND+k]<<"\n";
    }
  }
  In>>SLoad;                                                 //输入压力边界的个数
  if(SLoad>0)
  { SLoadNode=new int [SLoad][3];
    SLoadqo=new double [SLoad][3];
    for(i=0;i<SLoad;i++)
    { In>>SLoadNode[i][0]>>SLoadNode[i][1]>>SLoadNode[i][2];
      In>>SLoadqo[i][0]>>SLoadqo[i][1]>>SLoadqo[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)<<SLoadNode[i][j];
      for(j=0;j<3;j++)  Out<<setw(16)<<SLoadqo[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>>FixedNumber;                                               //输入位移约束个数
  if(FixedNumber<=0)
  { ofstream Error("file.err");
    Error<<"no input file ";Error.close();exit(1);
  }                                                   //未发现输入文件infile.dat
  FixedDOF=new int [FixedNumber][2];
  FixedData=new double [FixedNumber];

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

//完成所有数据的输入
 ///////////////////////////////////////////
 // 主体计算部分
////////////////////////////////////////////

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

  for(i=0;i<EleNumber;i++)
  { for(j=0;j<8;j++) 
      if(Ele[i][j]>0)
	  {   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,tt,xo,yo,ek,Ele[i],IK,nGauss);            //计算单元刚度矩阵
    Skke(ek,Je,16,zk,JD);                                 //组装单元刚度矩阵
    if(MLoad)
    CacuBF(Density,tt,ax,ay,xo,yo,F,Ele[i],IK,nGauss); }    //计算体积力载荷
                               //完成单元刚度矩阵的计算、组装以及体积力的计算
  
  }
  for(i=0;i<SLoad;i++)
  { for(j=0;j<3;j++) 
       if(SLoadNode[i][j])
	   {    xo[j]=X[SLoadNode[i][j]];
	        yo[j]=Y[SLoadNode[i][j]];
	   }                                                //返回当前单元边界坐标
    CacuSF(tt,SLoadqo[i],xo,yo,F,SLoadNode[i],IK,nGauss);  //计算表面力载荷
  }                                                     //完成表面分布载荷的处理
  for(i=0;i<FixedNumber;i++)                                         //处理边界条件
  Fixed((FixedDOF[i][0]-1)*ND+FixedDOF[i][1],FixedData[i],zk,F,JD,NodeNumber*ND);

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

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

  for(i=0;i<EleNumber;i++)
  {    for(j=0;j<8;j++)
          if(Ele[i][j]>=0)
		  {  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)
      Cass(E,u,xo,yo,uo,vo,Stress,Ele[i],IK);
     
  }
  for(i=1;i<NodeNumber;i++) for(j=1;j<5;j++)
  Stress[i][j]/=Stress[i][0];  // 应力绕节点平衡

  Out<<"\n                        节点应力         \n";
  Out<<"\n                  节点坐标\n";
  Out<<"  节点号     x向                   y向                  xy向                    z向\n";
  for(i=1;i<=NodeNumber;i++)
  Out<<setw(6)<<i<<setw(21)<<Stress[i][1]<<setw(21)<<Stress[i][2]<<setw(21)
  <<Stress[i][3]<<setw(21)<<Stress[i][4]<<"\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(SLoad>0) { delete [] SLoadNode; delete [] SLoadqo;}
  if(FixedDOF)  {delete [] FixedDOF;delete [] FixedData;}
  delete [] Stress;
  return ;
}