/*===========================*/ /* gccFRAME_EV.c */ /*===========================*/ #include #include #include #include #include #define NDmax 500 /*最大節点数*/ #define NEmax 500 /*最大要素数*/ #define MAmax 10 /*最大材料種別数*/ #define NTmax NDmax*3+1 /*最大自由度数*/ int main_part(char fnameR[50],char fnameW[50]); void DATAOUT_ISA(char fnameW[50],int NODT,int NELT, double x[NDmax+1],double y[NDmax+1],double dis[NTmax+1], double df[NTmax+1],double ftvec[NTmax+1],double reforce[NEmax+1][7]); void DATAOUT_EVA(char fnameW[50],int NODT,int NELT, double x[NDmax+1],double y[NDmax+1],double ev[NTmax+1],double vec[NTmax+1][NTmax+1]); void FRAMESM(int ne,int kakom[NEmax+1][3],double x[NDmax+1],double y[NDmax+1], double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double reforce[NEmax+1][7], double sme[7][7],double smg[7][7]); void ELMSM(int ne,double xx1,double yy1,double xx2,double yy2, double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double sme[7][7]); void CALST(int ne, int kakom[NEmax+1][3], double x[NDmax+1],double y[NDmax+1],double dis[NTmax+1], double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double fevec[7],double reforce[NEmax+1][7]); void HENMX(int ne,double xx1,double yy1,double xx2,double yy2,double hen[7][7]); void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][3],double sm[7][7],double tsm[NTmax+1][NTmax+1]); void ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][3],double fevec[7],double ftvec[NTmax+1]); int CHBAND10(int n1,int m1,double array[NTmax+1][NTmax+1]); void CHBAND11(int n1,int m1,double array[NTmax+1][NTmax+1],double vector[NTmax+1]); int IBAND(int mm,double array[NTmax+1][NTmax+1],double fact0); void BAND(int mm,int ib,double array[NTmax+1][NTmax+1]); void GJACOBI(int nnf,double a[NTmax+1][NTmax+1],double b[NTmax+1][NTmax+1],double ev[NTmax+1],double vec[NTmax+1][NTmax+1]); /*--------------------------------------------------------------------------------*/ int main(int argc, char* argv[]) { char fnameR[50]; /*入力ファイル名*/ char fnameW[50]; /*出力ファイル名*/ strcpy(fnameR,argv[1]); strcpy(fnameW,argv[2]); main_part(fnameR,fnameW); /* printf("リターンキーを押してください"); getchar(); */ return 0; } /*--------------------------------------------------------------------------------*/ int main_part(char fnameR[50],char fnameW[50]) { /*平面骨組固有値解析*/ FILE *fin,*fout; char dat[1000]=""; char sv[50]=""; int i,j,k,l,m,ne,n,ia,ja; char strcom[50]=""; /*コメント*/ int NODT; /*節点総数*/ int NELT; /*要素総数*/ int MATEL; /*材料種別数*/ int KOX; /*x方向変位既知節点数*/ int KOY; /*y方向変位既知節点数*/ int KOZ; /*回転変位既知節点数*/ int NF; /*載荷節点数*/ int nt; /*全自由度*/ int mm; /*変位既知節点処理語の自由度*/ int ib; /*バンド幅*/ double fact0=1.0e-30; /*0判定基準*/ clock_t t_sta,t_end; double cal_time; time_t jikoku; struct tm *lt; char *p; /*--------------------------------------------------------------------------------------*/ /*配列寸法宣言*/ /*--------------------------------------------------------------------------------------*/ int kakom[NEmax+1][3]; /*要素構成節点番号*/ double Em[NEmax+1]; /*要素弾性係数*/ double AA[NEmax+1]; /*要素断面積*/ double AI[NEmax+1]; /*要素断面二次モーメント*/ int matno[NEmax+1]; /*材料種別No*/ double wm[MAmax+1][4]; /*材料種別作業用*/ double x[NDmax+1]; /*初期節点x座標*/ double y[NDmax+1]; /*初期節点y座標*/ int nokx[NDmax+1]; /*x方向変位既知節点番号*/ int noky[NDmax+1]; /*y方向変位既知節点番号*/ int nokz[NDmax+1]; /*回転変位既知節点番号*/ int index[NTmax+1]; /*変位既知節点処理のためのインデックス*/ static double tsm[NTmax+1][NTmax+1]; /*線形全体剛性マトリックス*/ static double tsg[NTmax+1][NTmax+1]; /*非線形全体剛性マトリックス*/ static double vec[NTmax+1][NTmax+1]; /*固有ベクトル*/ double df[NTmax+1]; /*全体節点外力ベクトル増分(入力値)*/ double ftvec[NTmax+1]; /*全体節点力ベクトル*/ double dis[NTmax+1]; /*全体節点変位増分*/ double reac[NTmax+1]; /*全体節点変位増分*/ double ev[NTmax+1]; /*固有値*/ double sme[7][7]; /*要素剛性マトリックス*/ double smg[7][7]; /*要素剛性マトリックス*/ double fevec[7]; /*要素節点荷重ベクトル*/ double hen[7][7]; /*座標変換行列*/ double reforce[NEmax+1][7]; /*断面力*/ int jdg; /*特異行列判定*/ /*時間計測開始*/ t_sta=clock(); /*******************************************************/ /*データ入力*/ /*******************************************************/ fin=fopen(fnameR,"r"); /*--------------------------------------------------------------------------------------*/ /*書出用コメント入力*/ /*--------------------------------------------------------------------------------------*/ fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1); /*--------------------------------------------------------------------------------------*/ /*節点数,要素数,材料種別数][X方向変位既知節点数,Y方向変位既知節点数,回転変位既知節点数,載荷点数,初期弧長*/ /*--------------------------------------------------------------------------------------*/ fgets(dat,sizeof dat,fin); NODT =atoi(strtok(dat,",")); NELT =atoi(strtok(NULL,",")); MATEL =atoi(strtok(NULL,",")); KOX =atoi(strtok(NULL,",")); KOY =atoi(strtok(NULL,",")); KOZ =atoi(strtok(NULL,",")); NF =atoi(strtok(NULL,",")); nt=NODT*3; /*------------------------------------------------------------------------------------------*/ /*材料特性(Em,AA,AI)*/ /*------------------------------------------------------------------------------------------*/ for(i=1;i<=MATEL;i++){ fgets(dat,sizeof dat,fin); wm[i][1]=atof(strtok(dat,",")); wm[i][2]=atof(strtok(NULL,",")); wm[i][3]=atof(strtok(NULL,",")); } /*------------------------------------------------------------------------------------------*/ /*要素構成節点と材料種別(node-1,node-2,matno)*/ /*------------------------------------------------------------------------------------------*/ for(ne=1;ne<=NELT;ne++){ fgets(dat,sizeof dat,fin); kakom[ne][1]=atoi(strtok(dat,",")); kakom[ne][2]=atoi(strtok(NULL,",")); matno[ne] =atoi(strtok(NULL,",")); for(i=1;i<=MATEL;i++){ if(i==matno[ne]){ Em[ne] =wm[i][1]; AA[ne] =wm[i][2]; AI[ne] =wm[i][3]; } } } /*--------------------------------------------------------------------------------------*/ /*節点座標(x,y)*/ /*--------------------------------------------------------------------------------------*/ for(i=1;i<=NODT;i++){ fgets(dat,sizeof dat,fin); x[i] =atof(strtok(dat,",")); y[i] =atof(strtok(NULL,",")); } /*--------------------------------------------------------------------------------------*/ /*変位既知節点番号,既知節点変位量*/ /*--------------------------------------------------------------------------------------*/ if(00){mm=mm+1;index[mm]=i;}} /*******************************************************/ /*初期応力計算 {dis}=[K]^(-1){df}*/ /*******************************************************/ for(ne=1;ne<=NELT;ne++){for(i=1;i<=6;i++){reforce[ne][i]=0.0;}} for(i=1;i<=nt;i++){for(j=1;j<=nt;j++){tsm[i][j]=0.0;}} for(ne=1;ne<=NELT;ne++){ FRAMESM(ne,kakom,x,y,Em,AA,AI,reforce,sme,smg); ASMAT(ne,2,3,kakom,sme,tsm); } for(i=1;i<=nt;i++){reac[i]=0.0;dis[i]=0.0;ftvec[i]=0.0;} for(i=1;i<=mm;i++){ ia=index[i];reac[i]=df[ia]; for(j=1;j<=mm;j++){ja=index[j];tsm[i][j]=tsm[ia][ja];} } ib=IBAND(mm,tsm,fact0); BAND(mm,ib,tsm); jdg=CHBAND10(mm,ib,tsm);/*三角行列*/ if(jdg==1)goto epart; CHBAND11(mm,ib,tsm,reac);/*前進消去・後退代入*/ for(i=1;i<=mm;i++){ia=index[i];dis[ia]=reac[i];} /*内力計算*/ for(ne=1;ne<=NELT;ne++){ CALST(ne,kakom,x,y,dis,Em,AA,AI,fevec,reforce); ASVEC(ne,2,3,kakom,fevec,ftvec); } DATAOUT_ISA(fnameW,NODT,NELT,x,y,dis,df,ftvec,reforce); /*******************************************************/ /*固有値解析 [KL]{x}=ev[KG]{x}*/ /*******************************************************/ for(i=1;i<=nt;i++){for(j=1;j<=nt;j++){tsm[i][j]=0.0;tsg[i][j]=0.0;}} for(ne=1;ne<=NELT;ne++){ FRAMESM(ne,kakom,x,y,Em,AA,AI,reforce,sme,smg); ASMAT(ne,2,3,kakom,sme,tsm); ASMAT(ne,2,3,kakom,smg,tsg); } for(i=1;i<=nt;i++){ ev[i]=0.0; for(j=1;j<=nt;j++){vec[i][j]=0.0;} } for(i=1;i<=mm;i++){ ia=index[i]; for(j=1;j<=mm;j++){ ja=index[j]; tsm[i][j]=tsm[ia][ja]; tsg[i][j]=tsg[ia][ja]; } } GJACOBI(mm,tsm,tsg,ev,vec); for(i=1;i<=mm;i++){ printf("%g\n",ev[i]); } for(i=1;i<=nt;i++){ dis[i]=0.0; for(j=1;j<=nt;j++){tsm[i][j]=0.0;} } for(i=1;i<=mm;i++){ ia=index[i];dis[ia]=ev[i]; for(j=1;j<=mm;j++){ ja=index[j];tsm[ja][ia]=vec[j][i]; }} DATAOUT_EVA(fnameW,NODT,NELT,x,y,dis,tsm); epart: /*計算完了*/ t_end=clock(); cal_time=(double)(t_end-t_sta)/CLOCKS_PER_SEC; time(&jikoku); lt=localtime(&jikoku); p=asctime(lt); fout=fopen(fnameW,"a+"); fprintf(fout,"NODT=%d nt=%d mm=%d ib=%d\n",NODT,nt,mm,ib); fprintf(fout,"Calculation time=%.3f(sec)\n",cal_time); fprintf(fout,"%s\n",p); fclose(fout); printf("NODT=%d nt=%d mm=%d ib=%d\n",NODT,nt,mm,ib); printf("計算完了\n"); printf("計算時間:%.3f sec\n",cal_time); printf("%s\n",p); return 0; } /*--------------------------------------------------------------------------------*/ void DATAOUT_ISA(char fnameW[50],int NODT,int NELT, double x[NDmax+1],double y[NDmax+1],double dis[NTmax+1], double df[NTmax+1],double ftvec[NTmax+1],double reforce[NEmax+1][7]) { FILE *fout; int i,ne; char dat[256]=""; char sv[50]=""; fout=fopen(fnameW,"a+"); fprintf(fout,"***** Initial stress analysist *****\n"); fprintf(fout,"*Node value\n"); fprintf(fout,"node,x-coord,y-coord,x-dir,y-dir,z-dir,ff-x,ff-y,ff-z,ftvec-x,ftvec-y,ftvec-z\n"); for(i=1;i<=NODT;i++){ sprintf(sv,"%d",i) ;strcpy(dat,sv); sprintf(sv,"%e",x[i]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",y[i]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",dis[3*i-2]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",dis[3*i-1]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",dis[3*i]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",df[3*i-2]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",df[3*i-1]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",df[3*i]) ;strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",ftvec[3*i-2]);strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",ftvec[3*i-1]);strcat(dat,",");strcat(dat,sv); sprintf(sv,"%e",ftvec[3*i]) ;strcat(dat,",");strcat(dat,sv); fprintf(fout,"%s\n",dat); } fprintf(fout,"*Stress-resultant\n"); fprintf(fout,"element,Ni,Si,Mi,Nj,Sj,Mj\n"); for(ne=1;ne<=NELT;ne++){ sprintf(sv,"%d",ne);strcpy(dat,sv); for(i=1;i<=6;i++){ sprintf(sv,"%e",reforce[ne][i]);strcat(dat,",");strcat(dat,sv); } fprintf(fout,"%s\n",dat); } fclose(fout); } /*--------------------------------------------------------------------------------*/ void DATAOUT_EVA(char fnameW[50],int NODT,int NELT, double x[NDmax+1],double y[NDmax+1],double ev[NTmax+1],double vec[NTmax+1][NTmax+1]) { FILE *fout; int i,j,k,ne,nt; char dat[256]=""; char sv[50]=""; nt=NODT*3; fout=fopen(fnameW,"a+"); fprintf(fout,"***** Eigen value analysis *****\n"); i=0;k=0; do{ i=i+1; if(1e-30m1){ mm=m1; }else{ mm=n1-i+1; } for(j=1;j<=mm;j++){ z=array[i][j]; if(j!=m1){ if(m1-j+11){ array[i][j]=z/array[i][1]; }else{ if(z<0)return 1; array[i][1]=sqrt(z); } }else{ array[i][j]=z/array[i][1]; } } } return 0; } /*--------------------------------------------------------------------------------*/ void CHBAND11(int n1,int m1,double array[NTmax+1][NTmax+1],double vector[NTmax+1]) { /*前進消去・後退代入*/ int i,j,mm; double z; vector[1]=vector[1]/array[1][1]; for(i=2;i<=n1;i++){ z=vector[i]; if(i>m1){ mm=m1; }else{ mm=i; } for(j=2;j<=mm;j++){ z=z-array[i-j+1][j]*vector[i-j+1]; } vector[i]=z/array[i][1]; } vector[n1]=vector[n1]/array[n1][1]; for(i=n1-1;i>=1;i--){ z=vector[i]; for(j=2;j<=m1;j++){ if(i+j-1>n1)break; z=z-array[i][j]*vector[i+j-1]; } vector[i]=z/array[i][1]; } } /*--------------------------------------------------------------------------------*/ int IBAND(int mm,double array[NTmax+1][NTmax+1],double fact0) { int ib,i,j; /*-----------------------------------------------------*/ /*バンド幅計算 /*-----------------------------------------------------*/ ib=1; for(i=1;i<=mm-1;i++){ for(j=mm;j>=i+1;j--){ if(fact0fmk)fmk=fabs(a[i][j]); if(fabs(b[i][j])>fmg)fmg=fabs(b[i][j]); } } epsk=fmk*fact; epsg=fmg*fact; for(i=1;i<=nnf;i++){ for(j=1;j<=nnf;j++){ vec[i][j]=0.0; } vec[i][i]=1.0; } kmax=10*nnf*nnf; for(k=1;k<=kmax;k++){ ip=1; iq=2; if(fabs(a[1][1]*a[1][1]+a[2][2]*a[2][2])=amax){ im=i;amax=ev[i]; } } ww=ev[k];ev[k]=ev[im];ev[im]=ww; for(i=1;i<=nnf;i++){ ww=vec[i][k]; vec[i][k]=vec[i][im]; vec[i][im]=ww; } } } /*--------------------------------------------------------------------*/