/*===========================*/
/* gccFEM4nodASTS.c          */
/*===========================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 5000
#define NEmax 3000
#define MAmax 10
#define NTmax NDmax*2
#define pi M_PI

int main_part(char fnameR[50],char fnameW[50]);
void DMAT_AX4(double Eme,double pzr,double pot,double d[5][5]);
void BMAT_AX4(int ne,int kakom[NEmax+1][5],double z[NDmax+1],double r[NDmax+1],
    double bm[5][9],double *detJ,double a,double b);
void IPOINT4(int kk,double *a,double *b);
void SM_AX4(int ne,int kakom[NEmax+1][5],double z[NDmax+1],double r[NDmax+1],
    double Em[NEmax+1],double po[NEmax+1],double sm[9][9],int noten[NEmax+1][5]);
void CALST_AX4(int ne,int kakom[NEmax+1][5],double z[NDmax+1],double r[NDmax+1],
    double deltaT[NDmax+1],double dist[NTmax+1],
    double Em[NEmax+1],double po[NEmax+1],double alpha[NEmax+1],double ts[NEmax+1],
    double strsave[NEmax+1][5][5],double fevec[9],int noten[NEmax+1][5]);
void PST_AX4(double sigz,double sigr,double tauxy,double *ps1,double *ps2,double *ang);
void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][5],double sm[9][9],double tsm[NTmax+1][NTmax+1]);
void ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][5],double fevec[9],double ftvec[NTmax+1]);
void 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 a[NTmax+1][NTmax+1],double fact0);
void BAND(int mm,int ib,double a[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);

    return 0;
}
/*------------------------------------------------------------------------------*/

int main_part(char fnameR[50],char fnameW[50])
{
    /*四角形要素軸対称応力解析 */
    FILE *fin,*fout;
    char dat[256]="";
    char sv[20]="";
    int i,j,k,l,kk,ne,n,ia,ja,mn;
    char strcom[256]="";  /* 書出用コメント */
    int NODT;             /* 節点総数 */
    int NELT;             /* 要素総数 */
    int MATEL;            /* 材料種類数 */
    int KOZ;              /* z方向変位既知節点数 */
    int KOR;              /* r方向変位既知節点数 */
    int NF;               /* 載荷節点数 */
    int nt;               /* 全自由度（総節点数×２[１節点自由度]） */
    int nod=4;            /* 1要素節点数 */
    int nhen=2;           /* 1節点自由度数 */
    int mm;               /* 既知節点変位処理後自由度（逆行列の寸法） */
    int ib;               /* バンド幅 */
    double fact0=1.0E-30; /* 0判定 */
    double elarea1;       /* 要素半面積 */
    double elarea2;       /* 要素半面積 */
    double sigz,sigr,sigt,tauzr,ps1,ps2,ang;
    int IPR=0;            /* 0:全積分点出力,1:要素平均応力出力 */
    clock_t t_sta,t_end;
    double cal_time;
    time_t jikoku;
    struct tm *lt;
    char *p;
    /*--------------------------------------------------------------------------------------*/
    /*配列寸法宣言*/
    /*--------------------------------------------------------------------------------------*/
    int    kakom[NEmax+1][5];         /* 要素構成節点番号 */
    double Em[NEmax+1];               /* 要素弾性係数 */
    double po[NEmax+1];               /* 要素ポアソン比 */
    double alpha[NEmax+1];            /* 要素線膨張係数 */
    double ts[NEmax];                 /* 材料引張強度 */
    int    matno[NEmax+1];            /* 材料種別No */
    double wm[MAmax+1][5];            /* 材料物性作業用 */
    double strsave[NEmax+1][5][5];    /* 応力記憶用 */
    int    noten[NEmax+1][5];           /* No-tension積分点識別 */
    int    ntsum;                     /* No-tension積分点数 */
    int    nnn;                       /* 収束計算回数制御 */
    double dtest;                     /* 変位増分総和 */
    double ftest;                     /* 不平衡力総和 */
    int    icount;                    /* 規定以下変位増分自由度総和 */
    int    maxita=2000;               /* 収束計算回数上限 */
    double z[NDmax+1];                /* 節点z座標 */
    double r[NDmax+1];                /* 節点r座標 */
    double deltaT[NDmax+1];           /* 節点温度変化量[＋：温度上昇] */
    int    nokz[NDmax+1];             /* z方向変位既知節点番号 */
    int    nokr[NDmax+1];             /* r方向変位既知節点番号 */
    int    index[NTmax+1];          /* 既知節点変位処理のためのインデックス */
    static double tsm[NTmax+1][NTmax+1]; /* 全体剛性マトリックス */
    double f[NTmax+1];              /* 全体外力ベクトル */
    double ftvec[NTmax+1];          /* 全体荷重ベクトル */
    double dis[NTmax+1];            /* 全体節点変位増分 */
    double rdis[NTmax+1];            /* 全体節点変位増分 */
    double wdis[NTmax+1];            /* 全体節点変位増分 */
    double dist[NTmax+1];            /* 全体節点変位増分 */
    double reac[NTmax+1];           /* 作業用配列 */
    double sm[9][9];                  /* 要素剛性マトリックス */
    double fevec[9];                  /* 要素荷重ベクトル */

    /*時間計測開始*/
    t_sta=clock();
    /******************************************************/
    /*データ入力*/
    /******************************************************/
    fin=fopen(fnameR,"r");
    /*--------------------------------------------------------------------------------------*/
    /*書出用コメント入力*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1);
printf("%s\n",strcom);
    /*--------------------------------------------------------------------------------------*/
    /*応力状態[0;平面ひずみ]，節点数，要素数，材料種類数,X方向変位既知節点数，Y方向変位既知節点数，載荷点数*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);
    NODT  =atoi(strtok(dat,","));
    NELT  =atoi(strtok(NULL,","));
    MATEL =atoi(strtok(NULL,","));
    KOZ   =atoi(strtok(NULL,","));
    KOR   =atoi(strtok(NULL,","));
    NF    =atoi(strtok(NULL,","));
    IPR   =atoi(strtok(NULL,","));
    nt=NODT*2; /* 全自由度 */
printf("%d %d %d %d %d %d\n",NODT,NELT,MATEL,KOZ,KOR,NF);
    /*--------------------------------------------------------------------------------------*/
    /*材料特性（Em,po,alpha,ts）*/
    /*--------------------------------------------------------------------------------------*/
    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,","));
        wm[i][4]=atof(strtok(NULL,","));
    }
    /*--------------------------------------------------------------------------------------*/
    /*要素構成節点と材料種別（node-1,node-2,node-3,node-4,matno）*/
    /*--------------------------------------------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        fgets(dat,sizeof dat,fin);
        kakom[ne][1]=atoi(strtok(dat,","));
        kakom[ne][2]=atoi(strtok(NULL,","));
        kakom[ne][3]=atoi(strtok(NULL,","));
        kakom[ne][4]=atoi(strtok(NULL,","));
        matno[ne]=atoi(strtok(NULL,","));
        for(i=1;i<=MATEL;i++){
            if(i==matno[ne]){
                Em[ne]   =wm[i][1];
                po[ne]   =wm[i][2];
                alpha[ne]=wm[i][3];
                ts[ne]   =wm[i][4];
            }
        }
    }
    /*--------------------------------------------------------------------------------------*/
    /*節点座標（x,y），温度変化量*/
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=NODT;i++){
        fgets(dat,sizeof dat,fin);
        z[i]     =atof(strtok(dat,","));
        r[i]     =atof(strtok(NULL,","));
        deltaT[i]=atof(strtok(NULL,","));
    }
    /*--------------------------------------------------------------------------------------*/
    /*変位既知節点番号，既知節点変位量*/
    /*--------------------------------------------------------------------------------------*/
    if(0<KOZ){
        for(i=1;i<=KOZ;i++){
            fgets(dat,sizeof dat,fin);
            nokz[i]=atoi(strtok(dat,","));
            rdis[2*nokz[i]-1]=atof(strtok(NULL,","));
        }
    }
    if(0<KOR){
        for(i=1;i<=KOR;i++){
            fgets(dat,sizeof dat,fin);
            nokr[i]=atoi(strtok(dat,","));
            rdis[2*nokr[i]]=atof(strtok(NULL,","));
        }
    }
    /*--------------------------------------------------------------------------------------*/
    /* 節点番号,z方向荷重,r方向荷重 */
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=nt;i++){f[i]=0.0;}
    if(0<NF){
        for(i=1;i<=NF;i++){
            fgets(dat,sizeof dat,fin);
            n       =atoi(strtok(dat,","));
            f[2*n-1]=atof(strtok(NULL,","));
            f[2*n]  =atof(strtok(NULL,","));
        }
    }
    fclose(fin);
    /******************************************************/
    /* 入力確認出力 */
    /******************************************************/
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s\n",strcom);
    fprintf(fout,"NODT,NELT,MATEL,KOZ,KOR,NF,IPR\n");
    sprintf(sv,"%d",NODT) ;strcpy(dat,sv);
    sprintf(sv,"%d",NELT)   ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",MATEL)  ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOZ)    ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOR)    ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",NF)     ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",IPR)    ;strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
    fprintf(fout,"*node characteristics\n");
    fprintf(fout,"node,z,r,fz,fr,fix-z,fix-r,rdis-z,rdis-r,deltaT\n");
    for(i=1;i<=NODT;i++){
        k=0;l=0;
        if(0<KOZ){
            for(j=1;j<=KOZ;j++){
                if(i==nokz[j])k=1;
            }
        }
        if(0<KOR){
            for(j=1;j<=KOR;j++){
                if(i==nokr[j])l=1;
            }
        }
        sprintf(sv,"%d",i)           ;strcpy(dat,sv);
        sprintf(sv,"%e",z[i])        ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",r[i])        ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",f[2*i-1])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",f[2*i])      ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",k)           ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",l)           ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",rdis[2*i-1]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",rdis[2*i])   ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",deltaT[i])   ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fprintf(fout,"*element characteristics\n");
    fprintf(fout,"element,node-1,node-2,node-3,node-4,E,po,alpha,ts,matno\n");
    for(ne=1;ne<=NELT;ne++){
        sprintf(sv,"%d",ne)           ;strcpy(dat,sv);
        sprintf(sv,"%d",kakom[ne][1]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",kakom[ne][2]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",kakom[ne][3]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",kakom[ne][4]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",Em[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",po[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",alpha[ne])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ts[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",matno[ne])    ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fclose(fout);
    /*要素面積確認と異常終了通知*/
    for(ne=1;ne<=NELT;ne++){
        i=kakom[ne][1];j=kakom[ne][2];k=kakom[ne][3];
        elarea1=0.5*((z[k]-z[j])*r[i]+(z[i]-z[k])*r[j]+(z[j]-z[i])*r[k]);
        if(elarea1<fact0){
            /* 結果出力（Appendモード）*/
            fout=fopen(fnameW,"a+");
            fprintf(fout,"area of element is less or equal to 0\n");
            fclose(fout);
            printf("要素面積<=0：異常終了通知");
            fclose(fout);
            return 0;
        }
        i=kakom[ne][1];j=kakom[ne][3];k=kakom[ne][4];
        elarea2=0.5*((z[k]-z[j])*r[i]+(z[i]-z[k])*r[j]+(z[j]-z[i])*r[k]);
        if(elarea2<fact0){
            /* 結果出力（Appendモード）*/
            fout=fopen(fnameW,"a+");
            fprintf(fout,"area of element is less or equal to 0\n");
            fclose(fout);
            printf("要素面積<=0：異常終了通知");
            fclose(fout);
            return 0;
        }
    }
    /******************************************************/
    /*メイン処理計算*/
    /******************************************************/
    /*-----------------------------------------------------*/
    /*保存応力クリア*/
    /*-----------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        for(kk=1;kk<=4;kk++){
            for(i=1;i<=4;i++){strsave[ne][kk][i]=0.0;}
            noten[ne][kk]=0;
    }}
    /*-----------------------------------------------------*/
    /*全体剛性行列・全体荷重ベクトルクリア*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        ftvec[i]=0.0;
        wdis[i]=0.0;
        dist[i]=0.0;
        reac[i]=0.0;
        for(j=1;j<=nt;j++){
            tsm[i][j]=0.0;
        }
    }
    /*-----------------------------------------------------*/
    /*要素剛性行列・荷重ベクトル組立*/
    /*-----------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        /*剛性行列*/
        SM_AX4(ne,kakom,z,r,Em,po,sm,noten);
        ASMAT(ne,nod,nhen,kakom,sm,tsm);
    }
    /*-----------------------------------------------------*/
    /*荷重ベクトルセット*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        ftvec[i]=f[i];
    }
    /*-----------------------------------------------------*/
    /*変位既知節点自由度処理*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        index[i]=i;
    }
    if(0<KOZ){
        for(i=1;i<=KOZ;i++){n=2*nokz[i]-1;index[n]=0;}
    }
    if(0<KOR){
        for(i=1;i<=KOR;i++){n=2*nokr[i];index[n]=0;}
    }
    mm=0;
    for(i=1;i<=nt;i++){
        if(index[i]>0){
            mm=mm+1;index[mm]=i;
        }
    }
    for(i=1;i<=mm;i++){
        ia=index[i];
        ftvec[i]=ftvec[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);
    /*-----------------------------------------------------*/
    /*対角項確認と異常終了通知*/
    /*-----------------------------------------------------*/
    for(i=1;i<=mm;i++){
        if(fabs(tsm[i][1])<fact0){
            fout=fopen(fnameW,"a+");
            fprintf(fout,"diagonal value is less or equal to 0\n");
            fclose(fout);
            printf("対角項<=0：異常終了通知");
            fclose(fout);
            return 0;
        }
    }
    /*-----------------------------------------------------*/
    /*連立一次方程式（コレスキー法）*/
    /*-----------------------------------------------------*/
    CHBAND10(mm,ib,tsm);      /*三角行列*/
    /******************************************************/
    /*No-tension収束計算*/
    /******************************************************/
    nnn=0;
    do{
        nnn=nnn+1;
        /*-----------------------------------------------------*/
        /*連立一次方程式（コレスキー法）*/
        /*-----------------------------------------------------*/
        for(i=1;i<=mm;i++){
            reac[i]=ftvec[i];
        }
        CHBAND11(mm,ib,tsm,reac); /*前進消去・後退代入*/
        /*-----------------------------------------------------*/
        /*変位増分計算*/
        /*-----------------------------------------------------*/
        for(i=1;i<=nt;i++){dis[i]=0.0;dist[i]=0.0;}
        for(i=1;i<=mm;i++){ia=index[i];dis[ia]=reac[i];}
        for(i=1;i<=nt;i++){
            reac[i]=0.0;
            wdis[i]=wdis[i]+dis[i];
            dist[i]=wdis[i]+rdis[i];
        }
        /*-----------------------------------------------------*/
        /*応力計算・内力ベクトル組立*/
        /*-----------------------------------------------------*/
        for(ne=1;ne<=NELT;ne++){
            CALST_AX4(ne,kakom,z,r,deltaT,dist,Em,po,alpha,ts,strsave,fevec,noten);
            ASVEC(ne,nod,nhen,kakom,fevec,reac); /*内力ベクトル組立*/
        }
        /*-----------------------------------------------------*/
        /*不平衡力ベクトル算定*/
        /*-----------------------------------------------------*/
        for(i=1;i<=nt;i++){ftvec[i]=0.0;}
        for(i=1;i<=nt;i++){
            ftvec[i]=f[i]-reac[i];
        }
        for(i=1;i<=mm;i++){
            ia=index[i];
            ftvec[i]=ftvec[ia];
            dis[i]=dis[ia];
            dist[i]=dist[ia];
        }
        /*-----------------------------------------------------*/
        /*収束判定*/
        /*-----------------------------------------------------*/
        /*No-tension積分点の数*/
        ntsum = 0;
        for(ne=1;ne<=NELT;ne++){
            for(kk=1;kk<=4;kk++){
                ntsum = ntsum + fabs(noten[ne][kk]);
        }}
        /*変位増分絶対値・不平衡力絶対値総和*/
        icount=0;
        dtest=0.0;
        ftest=0.0;
        for(i=1;i<=mm;i++){
            dtest=dtest+fabs(dis[i]);   /*変位増分絶対値総和*/
            ftest=ftest+fabs(ftvec[i]); /*不平衡力絶対値総和*/
            if(fabs(dis[i]/dist[i])<1e-6)icount=icount+1;
        }
        printf("nnn=%d dtest=%e ftest=%e icount=%d\n",nnn,dtest,ftest,icount);

        /******************************************************/
        if(2<=nnn){
            if(ntsum==0)break;
            if(icount==mm)break;
        }
        /******************************************************/
    }while(nnn<=maxita);

    /* 出力用全変位調整 */
    for(i=1;i<=nt;i++){dis[i]=0.0;}
    for(i=1;i<=mm;i++){ia=index[i];dis[ia]=dist[i];}
    for(i=1;i<=nt;i++){dist[i]=dis[i]+rdis[i];}
    /******************************************************/
    /*結果出力*/
    /******************************************************/
    /*結果出力（Appendモード）*/
    fout=fopen(fnameW,"a+");
    fprintf(fout,"*displacements and forces\n");
    fprintf(fout,"node,coord-z,coord-r,dis-z,dis-r,reac-z,reac-r,ftvec-z,ftvec-r\n");
    /*節点変位*/
    for(i=1;i<=NODT;i++){
        sprintf(sv,"%d",i)          ;strcpy(dat,sv);
        sprintf(sv,"%e",z[i])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",r[i])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",dist[2*i-1]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",dist[2*i])  ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",reac[2*i-1]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",reac[2*i])  ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",f[2*i-1]-reac[2*i-1]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",f[2*i]  -reac[2*i])  ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    /*要素応力*/
    fprintf(fout,"*stresses\n");
    fprintf(fout,"element,kk,sig-z,sig-r,sig-t,tau-zr,ps1,ps2,ang,noten,matno\n");
    switch(IPR){
        case 0:
            for(ne=1;ne<=NELT;ne++){
                for(kk=1;kk<=4;kk++){
                    sigz=strsave[ne][kk][1];
                    sigr=strsave[ne][kk][2];
                    sigt=strsave[ne][kk][3];
                    tauzr=strsave[ne][kk][4];
                    PST_AX4(sigz,sigr,tauzr,&ps1,&ps2,&ang);
                    sprintf(sv,"%d",ne)           ;strcpy(dat,sv);
                    sprintf(sv,"%d",kk)           ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",sigz)         ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",sigr)         ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",sigt)         ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",tauzr)        ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",ps1)          ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",ps2)          ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%e",ang)          ;strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%d",noten[ne][kk]);strcat(dat,",");strcat(dat,sv);
                    sprintf(sv,"%d",matno[ne])    ;strcat(dat,",");strcat(dat,sv);
                    fprintf(fout,"%s\n",dat);
                }
            }
            break;
        case 1:
            for(ne=1;ne<=NELT;ne++){
                sigz=0.0;sigr=0.0;sigt=0.0;tauzr=0.0;ntsum=0;
                for(kk=1;kk<=4;kk++){
                    sigz =sigz +strsave[ne][kk][1];
                    sigr =sigr +strsave[ne][kk][2];
                    sigt =sigt +strsave[ne][kk][3];
                    tauzr=tauzr+strsave[ne][kk][4];
                    ntsum=ntsum+noten[ne][kk];
                }
                sigz =0.25*sigz;
                sigr =0.25*sigr;
                sigt =0.25*sigt;
                tauzr=0.25*tauzr;
                PST_AX4(sigz,sigr,tauzr,&ps1,&ps2,&ang);
                kk=0;
                sprintf(sv,"%d",ne)           ;strcpy(dat,sv);
                sprintf(sv,"%d",kk)           ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",sigz)         ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",sigr)         ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",sigt)         ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",tauzr)        ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",ps1)          ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",ps2)          ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%e",ang)          ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%d",ntsum)        ;strcat(dat,",");strcat(dat,sv);
                sprintf(sv,"%d",matno[ne])    ;strcat(dat,",");strcat(dat,sv);
                fprintf(fout,"%s\n",dat);
            }
            break;
    }
    /*計算完了*/
    t_end=clock();
    cal_time=(double)(t_end-t_sta)/CLOCKS_PER_SEC;
    time(&jikoku);
    lt=localtime(&jikoku);
    p=asctime(lt);
    fprintf(fout,"NODT=%d nt=%d mm=%d ib=%d\n",NODT,nt,mm,ib);
    fprintf(fout,"nnn=%d dtest=%e ftest=%e\n",nnn,dtest,ftest);
    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 DMAT_AX4(double Eme,double pzr,double pot,double d[5][5])
{
    /* 応力-ひずみ関係行列[De] */
    int i,j;

    for(i=1;i<=4;i++){
        for(j=1;j<=4;j++){
            d[i][j]=0.0;
        }
    }
    d[1][1]=1.0-pzr;
    d[1][2]=pzr;
    d[1][3]=pot;
    d[1][4]=0.0;
    d[2][1]=pzr;
    d[2][2]=1.0-pzr;
    d[2][3]=pot;
    d[2][4]=0.0;
    d[3][1]=pot;
    d[3][2]=pot;
    d[3][3]=1.0-pot;
    d[3][4]=0.0;
    d[4][1]=0.0;
    d[4][2]=0.0;
    d[4][3]=0.0;
    d[4][4]=0.5*(1.0-2.0*pzr);
    for(i=1;i<=4;i++){
        for(j=1;j<=4;j++){
            d[i][j]=d[i][j]*Eme/(1.0+pzr)/(1.0-2.0*pzr);
        }
    }
}
/*------------------------------------------------------------------------------*/

void BMAT_AX4(int ne,int kakom[NEmax+1][5],double z[NDmax+1],double r[NDmax+1],
    double bm[5][9],double *detJ,double a,double b)
{
    /* ひずみ-変位関係行列[B] */
    int i,j,k,l;
    double N1,N2,N3,N4,rr;
    double J11,J12,J21,J22;
    double dn1a,dn2a,dn3a,dn4a;
    double dn1b,dn2b,dn3b,dn4b;

    i=kakom[ne][1];j=kakom[ne][2];k=kakom[ne][3];l=kakom[ne][4];

    N1=0.25*(1.0-a)*(1.0-b);
    N2=0.25*(1.0+a)*(1.0-b);
    N3=0.25*(1.0+a)*(1.0+b);
    N4=0.25*(1.0-a)*(1.0+b);
    rr=N1*r[i]+N2*r[j]+N3*r[k]+N4*r[l];

    /* [dN/dadN/db] */
    dn1a=-0.25*(1.0-b);
    dn2a= 0.25*(1.0-b);
    dn3a= 0.25*(1.0+b);
    dn4a=-0.25*(1.0+b);
    dn1b=-0.25*(1.0-a);
    dn2b=-0.25*(1.0+a);
    dn3b= 0.25*(1.0+a);
    dn4b= 0.25*(1.0-a);

    /* Jacobiマトリックスとdet[J] */
    J11=dn1a*z[i]+dn2a*z[j]+dn3a*z[k]+dn4a*z[l];
    J12=dn1a*r[i]+dn2a*r[j]+dn3a*r[k]+dn4a*r[l];
    J21=dn1b*z[i]+dn2b*z[j]+dn3b*z[k]+dn4b*z[l];
    J22=dn1b*r[i]+dn2b*r[j]+dn3b*r[k]+dn4b*r[l];
    *detJ=J11*J22-J12*J21;

    /*[B]行列要素[dN/dxdN/dy] */
    for(i=1;i<=4;i++){
        for(j=1;j<=8;j++){
            bm[i][j]=0.0;
        }
    }
    bm[1][1]= J22*dn1a-J12*dn1b;
    bm[1][3]= J22*dn2a-J12*dn2b;
    bm[1][5]= J22*dn3a-J12*dn3b;
    bm[1][7]= J22*dn4a-J12*dn4b;

    bm[2][2]=-J21*dn1a+J11*dn1b;
    bm[2][4]=-J21*dn2a+J11*dn2b;
    bm[2][6]=-J21*dn3a+J11*dn3b;
    bm[2][8]=-J21*dn4a+J11*dn4b;

    bm[3][2]=N1/rr*(*detJ);
    bm[3][4]=N2/rr*(*detJ);
    bm[3][6]=N3/rr*(*detJ);
    bm[3][8]=N4/rr*(*detJ);

    bm[4][1]=-J21*dn1a+J11*dn1b;
    bm[4][2]= J22*dn1a-J12*dn1b;
    bm[4][3]=-J21*dn2a+J11*dn2b;
    bm[4][4]= J22*dn2a-J12*dn2b;
    bm[4][5]=-J21*dn3a+J11*dn3b;
    bm[4][6]= J22*dn3a-J12*dn3b;
    bm[4][7]=-J21*dn4a+J11*dn4b;
    bm[4][8]= J22*dn4a-J12*dn4b;

    for(i=1;i<=4;i++){
        for(j=1;j<=8;j++){
            bm[i][j]=bm[i][j]/(*detJ);
        }
    }

}
/*------------------------------------------------------------------------------*/

void IPOINT4(int kk,double *a,double *b)
{
    /* Gauss積分点 */
    *a=0.0;*b=0.0;
    switch(kk){
        case 1:*a=-0.5773502692;*b=-0.5773502692;break;
        case 2:*a= 0.5773502692;*b=-0.5773502692;break;
        case 3:*a= 0.5773502692;*b= 0.5773502692;break;
        case 4:*a=-0.5773502692;*b= 0.5773502692;break;
    }
}
/*------------------------------------------------------------------------------*/

void SM_AX4(int ne,int kakom[NEmax+1][5],double z[NDmax+1],double r[NDmax+1],
    double Em[NEmax+1],double po[NEmax+1],double sm[9][9],int noten[NEmax+1][5])
{
    /* 要素剛性行列 */
    int kk,i,j,k;
    int i1,i2,i3,i4;
    double d[5][5];
    double bm[5][9];
    double str[5][9];
    double detJ,a,b;
    double N1,N2,N3,N4,rr;

    /* 配列クリア */
    for(i=1;i<=8;i++){
        for(j=1;j<=8;j++){
            sm[i][j]=0.0;
        }
    }
    /* 剛性行列[sm]=[B]T[D][B]*r*detJ */
    DMAT_AX4(Em[ne],po[ne],po[ne],d);

    i1=kakom[ne][1];i2=kakom[ne][2];i3=kakom[ne][3];i4=kakom[ne][4];

    for(kk=1;kk<=4;kk++){
        IPOINT4(kk,&a,&b);
        BMAT_AX4(ne,kakom,z,r,bm,&detJ,a,b);
        for(i=1;i<=4;i++){
            for(j=1;j<=8;j++){
                str[i][j]=0.0;
                for(k=1;k<=4;k++){
                    str[i][j]=str[i][j]+d[i][k]*bm[k][j];
                }
            }
        }
        N1=0.25*(1.0-a)*(1.0-b);
        N2=0.25*(1.0+a)*(1.0-b);
        N3=0.25*(1.0+a)*(1.0+b);
        N4=0.25*(1.0-a)*(1.0+b);
        rr=N1*r[i1]+N2*r[i2]+N3*r[i3]+N4*r[i4];
        for(i=1;i<=8;i++){
            for(j=1;j<=8;j++){
                for(k=1;k<=4;k++){
                    sm[i][j]=sm[i][j]+bm[k][i]*str[k][j]*rr*detJ;
                }
            }
        }
    }

}
/*------------------------------------------------------------------------------*/

void CALST_AX4(int ne,int kakom[NEmax+1][5],double z[NDmax+1],double r[NDmax+1],
    double deltaT[NDmax+1],double dist[NTmax+1],
    double Em[NEmax+1],double po[NEmax+1],double alpha[NEmax+1],double ts[NEmax+1],
    double strsave[NEmax+1][5][5],double fevec[9],int noten[NEmax+1][5])
{
    /* 応力および内力ベクトル計算 */
    int kk,i,j,k,ia;
    int i1,i2,i3,i4;
    double detJ,a,b;
    double d[5][5];
    double bm[5][9];
    double epsilon[5];
    double stress[5];
    double wd[9];
    double N1,N2,N3,N4,rr;
    double tem,eps0[5],pzr,pot;
    double sigz,sigr,sigt,tauzr,ps1,ps2,ang,phi;

    i1=kakom[ne][1];i2=kakom[ne][2];i3=kakom[ne][3];i4=kakom[ne][4];

    /* 要素座標節点変位 */
    for(i=1;i<=4;i++){
        ia=kakom[ne][i];
        wd[2*i-1]=dist[2*ia-1];
        wd[2*i]=dist[2*ia];
    }

    /* 要素応力*/
    for(kk=1;kk<=4;kk++){
        IPOINT4(kk,&a,&b);
        BMAT_AX4(ne,kakom,z,r,bm,&detJ,a,b);
        /*ひずみ計算 {epsilon}=[B]{u} */
        for(i=1;i<=4;i++){
            epsilon[i]=0.0;
            for(j=1;j<=8;j++){
                epsilon[i]=epsilon[i]+bm[i][j]*wd[j];
            }
        }
        /*積分点温度ひずみ*/
        N1=0.25*(1.0-a)*(1.0-b);
        N2=0.25*(1.0+a)*(1.0-b);
        N3=0.25*(1.0+a)*(1.0+b);
        N4=0.25*(1.0-a)*(1.0+b);
        tem=N1*deltaT[i1]+N2*deltaT[i2]+N3*deltaT[i3]+N4*deltaT[i4];
        eps0[1]=alpha[ne]*tem;
        eps0[2]=alpha[ne]*tem;
        eps0[3]=alpha[ne]*tem;
        eps0[4]=0.0;
        switch(noten[ne][kk]){
            case 0:pzr=po[ne];pot=po[ne];break;
            case 1:pzr=0.0   ;pot=po[ne];break;
            case 2:pzr=0.0   ;pot=po[ne];break;
            case 3:pzr=po[ne];pot=0.0   ;break;
            case 4:pzr=0.0   ;pot=0.0   ;break;
            case 5:pzr=0.0   ;pot=0.0   ;break;
        }
        DMAT_AX4(Em[ne],pzr,pot,d);
        /*応力計算 {stress}=[D]{epsilon} */
        for(i=1;i<=4;i++){
            stress[i]=0.0;
            for(j=1;j<=4;j++){
                stress[i]=stress[i]+d[i][j]*(epsilon[j]-eps0[j]);
            }
        }

        /*現状応力保存*/
        sigz =stress[1];
        sigr =stress[2];
        sigt =stress[3];
        tauzr=stress[4];
        PST_AX4(sigz,sigr,tauzr,&ps1,&ps2,&ang);
        /*No-tension判定*/
        noten[ne][kk]=0;
        if(ts[ne]<ps1&&ps2<ts[ne]&&sigt<ts[ne])noten[ne][kk]=1;
        if(ts[ne]<ps1&&ts[ne]<ps2&&sigt<ts[ne])noten[ne][kk]=2;
        if(ps1<ts[ne]&&ps2<ts[ne]&&ts[ne]<sigt)noten[ne][kk]=3;
        if(ts[ne]<ps1&&ps2<ts[ne]&&ts[ne]<sigt)noten[ne][kk]=4;
        if(ts[ne]<ps1&&ts[ne]<ps2&&ts[ne]<sigt)noten[ne][kk]=5;
        /*引張成分除去*/
        phi=-ang/180.0*pi;
        if(noten[ne][kk]==1||noten[ne][kk]==2||noten[ne][kk]==4||noten[ne][kk]==5){
            sigz =sigz -0.5*ps1*(1.0+cos(2.0*phi));
            sigr =sigr -0.5*ps1*(1.0-cos(2.0*phi));
            tauzr=tauzr+0.5*ps1*sin(2.0*phi);
        }
        if(noten[ne][kk]==2||noten[ne][kk]==5){
            sigz =sigz -0.5*ps2*(1.0-cos(2.0*phi));
            sigr =sigr -0.5*ps2*(1.0+cos(2.0*phi));
            tauzr=tauzr-0.5*ps2*sin(2.0*phi);
        }
        if(noten[ne][kk]==3||noten[ne][kk]==4||noten[ne][kk]==5){
            sigt=0.0;
        }
        /*作業用応力保存*/
        strsave[ne][kk][1]=sigz;
        strsave[ne][kk][2]=sigr;
        strsave[ne][kk][3]=sigt;
        strsave[ne][kk][4]=tauzr;
    }

    /* 内力ベクトル */
    for(i=1;i<=8;i++){fevec[i]=0.0;}
    for(kk=1;kk<=4;kk++){
        IPOINT4(kk,&a,&b);
        BMAT_AX4(ne,kakom,z,r,bm,&detJ,a,b);
        N1=0.25*(1.0-a)*(1.0-b);
        N2=0.25*(1.0+a)*(1.0-b);
        N3=0.25*(1.0+a)*(1.0+b);
        N4=0.25*(1.0-a)*(1.0+b);
        rr=N1*r[i1]+N2*r[i2]+N3*r[i3]+N4*r[i4];
        for(i=1;i<=8;i++){
            for(j=1;j<=4;j++){
                fevec[i]=fevec[i]+bm[j][i]*strsave[ne][kk][j]*rr*detJ;
            }
        }
    }

}
/*------------------------------------------------------------------------------*/
void PST_AX4(double sigz,double sigr,double tauzr,double *ps1,double *ps2,double *ang)
{
    /*主応力計算*/
    *ps1=0.5*(sigz+sigr)+sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr);
    *ps2=0.5*(sigz+sigr)-sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr);
    *ang=0.0;
    if(sigz==sigr){
        if(tauzr>0)    *ang=45.0;
        if(tauzr<0)    *ang=135.0;
        if(tauzr==0.0) *ang=0.0;
    }else{
        *ang=180.0/pi*(0.5*atan(2.0*tauzr/(sigz-sigr)));
        if(sigz>sigr&&tauzr>=0) *ang=*ang+0.0;
        if(sigz>sigr&&tauzr<0)  *ang=*ang+180.0;
        if(sigz<sigr)           *ang=*ang+90.0;
    }

}
/*------------------------------------------------------------------------------*/

void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][5],double sm[9][9],double tsm[NTmax+1][NTmax+1])
{
    /*剛性行列組立[nod：1要素節点数，nhen：1節点の自由度]*/
    int i,j,k,l,ki,kj,ks,js,kik,kjl,ksk,jsl;

    for(i=1;i<=nod;i++){
        for(j=1;j<=nod;j++){
            ki=(kakom[ne][i]-1)*nhen;
            kj=(kakom[ne][j]-1)*nhen;
            ks=(i-1)*nhen;
            js=(j-1)*nhen;
            for(k=1;k<=nhen;k++){
                for(l=1;l<=nhen;l++){
                    kik=ki+k;
                    kjl=kj+l;
                    ksk=ks+k;
                    jsl=js+l;
                    tsm[kik][kjl]=tsm[kik][kjl]+sm[ksk][jsl];
                }
            }
        }
    }

}
/*------------------------------------------------------------------------------*/

void ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][5],double fevec[9],double ftvec[NTmax+1])
{
    /*荷重ベクトル組立[nod：1要素節点数，nhen：1節点の自由度]*/
    int i,k,ki,ks,kik,ksk;

    for(i=1;i<=nod;i++){
        ki=(kakom[ne][i]-1)*nhen;
        ks=(i-1)*nhen;
        for(k=1;k<=nhen;k++){
            kik=ki+k;
            ksk=ks+k;
            ftvec[kik]=ftvec[kik]+fevec[ksk];
        }
    }

}
/*------------------------------------------------------------------------------*/

void CHBAND10(int n1,int m1,double array[NTmax+1][NTmax+1])
{
    /*三角行列*/
    int i,j,k,mm,mn;
    double z;

    array[1][1]=sqrt(array[1][1]);
    for(j=2;j<=m1;j++){
        array[1][j]=array[1][j]/array[1][1];
    }
    for(i=2;i<=n1;i++){
        if(n1-i+1>m1){
            mm=m1;
        }else{
            mm=n1-i+1;
        }
        for(j=1;j<=mm;j++){
            z=array[i][j];
            if(j!=m1){
                if(m1-j+1<i){
                    mn=m1-j+1;
                }else{
                    mn=i;
                }
                for(k=2;k<=mn;k++){
                    z=z-array[i-k+1][k]*array[i-k+1][j+k-1];
                }
                if(j>1){
                    array[i][j]=z/array[i][1];
                }else{
                    array[i][1]=sqrt(z);
                }
            }else{
                array[i][j]=z/array[i][1];
            }
        }
    }
}
/*--------------------------------------------------------------------------------*/

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(fact0<fabs(array[i][j]))break;
        }
        if(ib<=j-i+1)ib=j-i+1;
    }
    if(mm<ib)ib=mm;
    return ib;
}
/*--------------------------------------------------------------------------------*/

void BAND(int mm,int ib,double array[NTmax+1][NTmax+1])
{
    int i,j,k,mn;
    double work[mm+1];
    /*-----------------------------------------------------*/
    /*フルマトリックスのバンド化記憶
    /*-----------------------------------------------------*/
    for(i=2;i<=mm;i++){
        if(i+ib-1<mm){
            mn=ib;
        }else{
            mn=mm-i+1;
        }
        for(j=1;j<=ib;j++){
            work[j]=0.0;
        }
        for(k=1;k<=mn;k++){
            work[k]=array[i][i+k-1];
        }
        for(j=1;j<=ib;j++){
            array[i][j]=work[j];
        }
    }
}
/*--------------------------------------------------------------------------------*/