interface_c_and_ANSYS.rar

#define ANSYS #include"string.h" #include"stdio.h" #include"conio.h" #include"stdlib.h" #include"process.h" #include"malloc.h" #include"math.h" #define PI 3.1415926 float NODE_X[500],NODE_Y[500],INFO_LOADLINE[30][5],INFO_LOADNODE[30][3],WK[4],*stiff; int ELE_NUM[1000][3],BC_INFO[50][3],NODE_SUM,ELE_SUM,BC_LINE,BC_NODE,LOAD_LINE,LOAD_NODE,PTYPE; int bc_num,load_num,NJ2,SH[8],LL[6],ID[1000]; float DISPLACE[1000],D[3][3],EK[6][6],S[3][6]; FILE *fp,*fp1; char Input(char *); //读入相关参数 float Asek(int n,int iask); //求解单元刚度矩阵 void Astk(),Astp(),Cons(),Oned(),Solve(); void Result(); //输出相关结果 void main() {char instr[30]; //存放输入和输出的文件名 int i; printf("\nplease input the datafile name:"); scanf("%s",instr); //输入原始数据文件名 i=Input(instr); //读入数据，返回信息变量 if(i==1){printf("Data file input OK!Any key continue\n");getch();} if(i==2){printf("Data file open error!\n");exit(0);} Astk(); //组装总刚矩阵 Astp(); //组装总载矩阵 Cons(); //约束处理 Solve(); //求解方程组 printf("\nCalculating finish.Anykey to continue...\n"); getch(); Result(); //给出结果并输出 printf("\nAnykey to QUIT!\n"); getch(); } char Input(char *str) {int i,j,nd,k1,k2,j1,n,sum,ss=0,sp=0; //sp标识当前已读入的信息数据个数个数 float x1,y1,x2,y2,nr,dx1,dy1,dk,dx,dyx,q1,q2,aL; float a1,a2,a3,a4,w2[20]; float *w1; //指向读入的信息数据的临时存储空间 if((fp=fopen(str,"r+"))==NULL)return 2; //打开数据文件 if(ferror(fp))return 2; fscanf(fp,"%d %d %d %d %d %d %d",&PTYPE,&NODE_SUM,&ELE_SUM,&BC_LINE,&BC_NODE,&LOAD_LINE,&LOAD_NODE); //读入控制数据，PTYPE为1表示平面应力问题，为2表示平面应变问题 #ifdef ANSYS //条件编译，若定义了“ANSYS”则执行下面的语句 sum=(3*NODE_SUM)+(14*ELE_SUM)+(5*BC_LINE)+(3*BC_NODE)+(7*LOAD_LINE)+(3*LOAD_NODE)+4; //计算需要读入的信息数据个数 if(!(w1=(float *)malloc(sum*sizeof(float)))) //开辟临时空间给信息数据，首地址为w1 {printf("out of memory.\n"); exit(0); } //若空间不够则报警 if((fp1=fopen("node_ansys.in","r+"))==NULL)return 2; for(ss=0;ss<3*NODE_SUM;ss++) fscanf(fp1,"%f",w1+ss); //从w1地址开始存入节点信息数据 fclose(fp1); if((fp1=fopen("element_ansys.in","r+"))==NULL)return 2; for(ss=3*NODE_SUM;ss<(3*NODE_SUM+14*ELE_SUM);ss++) fscanf(fp1,"%f",w1+ss); //继续存入单元信息数据 fclose(fp1); for(ss=(3*NODE_SUM+14*ELE_SUM);ss<sum;ss++) fscanf(fp,"%f",w1+ss); //继续存入其他信息数据 #else //条件编译，若没有定义“ANSYS”则执行下面的语句 sum=(3*NODE_SUM)+(4*ELE_SUM)+(5*BC_LINE)+(3*BC_NODE)+(7*LOAD_LINE)+(3*LOAD_NODE)+4; if(!(w1=(float *)malloc(sum*sizeof(float)))) {printf("out of memory.\n"); exit(0); } if((fp1=fopen("node.in","r+"))==NULL)return 2; for(ss=0;ss<3*NODE_SUM;ss++) fscanf(fp1,"%f",w1+ss); fclose(fp1); if((fp1=fopen("element.in","r+"))==NULL)return 2; for(ss=3*NODE_SUM;ss<(3*NODE_SUM+4*ELE_SUM);ss++) fscanf(fp1,"%f",w1+ss); fclose(fp1); for(ss=(3*NODE_SUM+4*ELE_SUM);ss<sum;ss++) fscanf(fp,"%f",w1+ss); #endif //条件编译结束 fclose(fp); if(NODE_SUM>0) //NODE_SUM表示节点总数 {for(i=1;i<=NODE_SUM;i++) {nd=*(w1+sp+3*i-2-1); //nd读入节点编号 NODE_X[nd-1]=*(w1+sp+3*i-1-1); //NODE_X[nd-1]存储该节点横坐标 NODE_Y[nd-1]=*(w1+sp+3*i-0-1); } //NODE_Y[nd-1]存储该节点纵坐标 sp=sp+3*NODE_SUM;} if(ELE_SUM>0) //ELE_SUM表示单元总数 #ifdef ANSYS //条件编译，若定义了“ANSYS”则执行下面的语句 {for(i=1;i<=ELE_SUM;i++) {nd=*(w1+sp+14*i-1); //nd读入单元编号 ELE_NUM[nd-1][0]=*(w1+sp+14*i-13-1); ELE_NUM[nd-1][1]=*(w1+sp+14*i-12-1); ELE_NUM[nd-1][2]=*(w1+sp+14*i-11-1); } //ELE_NUM[nd-1]存储该单元三节点编码 sp=sp+14*ELE_SUM;} #else //条件编译，若没有定义“ANSYS”则执行下面的语句 {for(i=1;i<=ELE_SUM;i++) {nd=*(w1+sp+4*i-3-1); //nd读入单元编号 ELE_NUM[nd-1][0]=*(w1+sp+4*i-2-1); ELE_NUM[nd-1][1]=*(w1+sp+4*i-1-1); ELE_NUM[nd-1][2]=*(w1+sp+4*i-1); } //ELE_NUM[nd-1]存储该单元三节点编码 sp=sp+4*ELE_SUM;} #endif //条件编译结束 bc_num=0; if(BC_LINE>0) //BC_LINE表示有约束的直线边数 {for(i=1;i<=BC_LINE;i++) {x1=*(w1+sp+5*i-5);y1=*(w1+sp+5*i-4); x2=*(w1+sp+5*i-3);y2=*(w1+sp+5*i-2); //存储该边的起终点坐标 nr=*(w1+sp+5*i-1); //nr读入该边的约束信息，-1则x方向有约束，0则x,y方向均有约束，1则y方向有约束 k1=0;k2=0; if(nr<0.1) k1=1; //x方向有约束则令k1为1 if(nr>-0.1) k2=1; //y方向有约束则令k2为1 dx=x2-x1; //下面分斜率是否存在两种情况分别寻找在直线上的节点，并把边的约束信息转化到节点上 if(fabs(dx)<0.0001) //斜率不存在 {for(j=1;j<=NODE_SUM;j++) {dx1=NODE_X[j-1]-x1; if(fabs(dx1)<0.0001) //顺次判断各个节点是否在该直线上 {bc_num++; //存储数组元素加1 BC_INFO[bc_num-1][0]=j; BC_INFO[bc_num-1][1]=k1; BC_INFO[bc_num-1][2]=k2; } //BC_INFO[]存储在该直线上的节点的约束信息，依次为节点编号和x,y方向的约束情况 } } else //斜率存在 {dk=(y2-y1)/dx; //求斜率 for(j=1;j<=NODE_SUM;j++) {dy1=NODE_Y[j-1]-y1-dk*(NODE_X[j-1]-x1); if(fabs(dy1)<0.0001) //顺次判断各个节点是否在该直线上 {bc_num++; BC_INFO[bc_num-1][0]=j; BC_INFO[bc_num-1][1]=k1; BC_INFO[bc_num-1][2]=k2;} } } } sp=sp+5*BC_LINE;} if(BC_NODE>0) //BC_NODE表示有约束的节点数 {for(i=1;i<=BC_NODE;i++) {bc_num++; BC_INFO[bc_num-1][0]=*(w1+sp+3*i-2-1); BC_INFO[bc_num-1][1]=*(w1+sp+3*i-1-1); BC_INFO[bc_num-1][2]=*(w1+sp+3*i-1);} //BC_INFO[]存储节点约束信息，依次为节点编号和x,y方向的约束情况 sp=sp+3*BC_NODE;} load_num=0; if(LOAD_LINE>0) //LOAD_LINE表示有分布载荷的直线边数 {for(i=1;i<=LOAD_LINE;i++) {j1=0; x1=*(w1+sp+7*i-6-1);y1=*(w1+sp+7*i-5-1); x2=*(w1+sp+7*i-4-1);y2=*(w1+sp+7*i-3-1); //读入该边的约束始终点坐标 q1=*(w1+sp+7*i-2-1);q2=*(w1+sp+7*i-1-1); //读入始终点的载荷集度 aL=*(w1+sp+7*i-1); //读入从x轴按右手法则转至该载荷外作用线的角度 dx=x2-x1; if(fabs(dx)<0.0001) //斜率不存在 {for(j=1;j<=NODE_SUM;j++) {dx1=NODE_X[j-1]-x1; if(fabs(dx1)<0.0001) //顺次判断各节点是否在该直线上 {j1++; w2[j1-1]=j; } //w2[]存储在该直线上的节点编号 } } else //斜率存在 {dk=(y2-y1)/dx; //求斜率 for(j=1;j<=NODE_SUM;j++) {dy1=NODE_Y[j-1]-y1-dk*(NODE_X[j-1]-x1); if(fabs(dy1)<0.0001) //顺次判断各节点是否在该直线上 {j1++; //j1最后得到在该直线上的节点数 w2[j1-1]=j;} //w2[]存储在该直线上的节点编号 } } a1=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)); a2=(q2-q1)/a1; //求出载荷梯度 //下面一个循环将对直线边的分布载荷转化到各个单元边上 for(n=1;n<=j1-1;n++) {load_num++; k1=w2[n-1];k2=w2[n]; //单元边始终节点编码 a3=sqrt((NODE_X[k1-1]-x1)*(NODE_X[k1-1]-x1)+(NODE_Y[k1-1]-y1)*(NODE_Y[k1-1]-y1)); a4=sqrt((NODE_X[k2-1]-x1)*(NODE_X[k2-1]-x1)+(NODE_Y[k2-1]-y1)*(NODE_Y[k2-1]-y1)); //a3,a4分别保存始终点到约束起点的距离 INFO_LOADLINE[load_num-1][0]=k1; INFO_LOADLINE[load_num-1][1]=k2; //第load_num条单元边的始终节点编码 INFO_LOADLINE[load_num-1][2]=q1+a3*a2; INFO_LOADLINE[load_num-1][3]=q1+a4*a2; //第load_num条单元边的始终载荷 INFO_LOADLINE[load_num-1][4]=aL;} //第load_num条单元边的载荷方向 } sp=sp+7*LOAD_LINE;} if(LOAD_NODE>0) //LOAD_NODE表示集中力作用节点数 {for(i=1;i<=LOAD_NODE;i++) {INFO_LOADNODE[i-1][0]=*(w1+sp+3*i-2-1); INFO_LOADNODE[i-1][1]=*(w1+sp+3*i-1-1); INFO_LOADNODE[i-1][2]=*(w1+sp+3*i-1); } //INFO_LOADNODE[]存储集中力作用点信息，依次为节点坐标，力的大小，作用方向 sp=sp+3*LOAD_NODE;} for(i=0;i<4;i++)WK[i]=*(w1+sp+i); //WK[]依次读入材料的弹性模量,泊松比,容重和板厚 free(w1); //释放存储信息数据的临时空间 return 1;} float Asek(int n,int iask) //计算单元刚度矩阵，n为单元编号 {float bi,ci,cm,bm,cj,bj,b[3][6],bb[6][3]; int i,j,m,k1=0,k2=0,k3,k4; float th,ae,a2,at; float x1,x2,x3,y1,y2,y3; i=ELE_NUM[n-1][0];j=ELE_NUM[n-1][1];m=ELE_NUM[n-1][2]; //节点单元码i,j,m cm=NODE_X[j-1]-NODE_X[i-1];bm=NODE_Y[i-1]-NODE_Y[j-