/*===========================*/
/* gccFRAME.c                */
/*===========================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 500        /*最大節点数*/
#define NEmax 500        /*最大要素数*/
#define MAmax 200        /*最大材料種別数*/
#define Nnod   2         /*1要素節点数*/
#define Nfre   3         /*1節点自由度*/
#define Nfel  Nnod*Nfre  /*1要素自由度*/
#define NTmax NDmax*Nfre /*最大自由度数*/
int main_part(char fnameR[],char fnameW[]);
void SM_FRM(int ne,int kakom[NEmax+1][Nnod+1],double x[NDmax+1],double y[NDmax+1],
    double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double sm[Nfel+1][Nfel+1],double hen[Nfel+1][Nfel+1]);
void MV_FRM(int ne,int kakom[NEmax+1][Nnod+1],double x[NDmax+1],double y[NDmax+1],
    double AA[NEmax+1],double gamma[NEmax+1],double gkh[NEmax+1],double gkv[NEmax+1],
    double fevec[Nfel+1],double hen[Nfel+1][Nfel+1]);
void TV_FRM(int ne,int kakom[NEmax+1][Nnod+1],double deltaT[NDmax+1],
    double Em[NEmax+1],double AA[NEmax+1],double alpha[NEmax+1],double fevec[Nfel+1],double hen[Nfel+1][Nfel+1]);
void CALST_FRM(int ne, int kakom[NEmax+1][Nnod+1], double dis[NTmax+1], double sm[Nfel+1][Nfel+1], double hen[Nfel+1][Nfel+1], double reac[NTmax+1],
    double tsm[NTmax+1][NTmax+1],double Em[NEmax+1],double AA[NEmax+1],double alpha[NEmax+1],double deltaT[NDmax+1]);
void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][Nnod+1],double sm[Nfel+1][Nfel+1],double tsm[NTmax+1][NTmax+1]);
void ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][Nnod+1],double fevec[Nfel+1],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 array[NTmax+1][NTmax+1],double fact0);
void BAND(int mm,int ib,double array[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[],char fnameW[])
{
    /*平面骨組構造解析*/
    FILE *fin,*fout;
    char dat[256]="";
    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 nod;            /*1要素節点数*/
    int nhen;           /*1節点自由度*/
    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][Nnod+1];/*要素構成節点番号*/
    double Em[NEmax+1];           /*要素弾性係数*/
    double AA[NEmax+1];           /*要素断面積*/
    double AI[NEmax+1];           /*要素断面二次モーメント*/
    double gamma[NEmax+1];        /*要素単位体積重量*/
    double gkh[NEmax+1];          /*要素水平加速度（重力加速度比）*/
    double gkv[NEmax+1];          /*要素鉛直加速度（重力加速度比）*/
    double alpha[NEmax+1];        /*要素線膨張係数*/
    int    matno[NEmax+1];        /*材料種別No*/
    double wm[MAmax+1][8];        /*材料種別作業用*/
    double x[NDmax+1];            /*節点x座標*/
    double y[NDmax+1];            /*節点y座標*/
    double deltaT[NDmax+1];       /*節点温度変化量*/
    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]; /*全体剛性マトリックス*/
    double f[NTmax+1];            /*全体節点外力ベクトル（入力値）*/
    double fmass[NTmax+1];        /*全体節点慣性力ベクトル（計算値）*/
    double ftemp[NTmax+1];        /*全体節点温度荷重相当ベクトル（計算値）*/
    double ftvec[NTmax+1];        /*全体節点力ベクトル*/
    double fdis[NTmax+1];         /*既知節点変位による荷重成分格納*/
    double rdis[NTmax+1];         /*既知節点変位量*/
    double dis[NTmax+1];          /*全体節点変位*/
    double reac[NTmax+1];         /*要素節点反力*/
    double sm[Nfel+1][Nfel+1];    /*要素剛性マトリックス*/
    double fevec[Nfel+1];         /*要素節点荷重ベクトル*/
    double hen[Nfel+1][Nfel+1];   /*座標変換行列*/

    /*時間計測開始*/
    t_sta=clock();
    /*******************************************************/
    /*データ入力*/
    /*******************************************************/
    fin=fopen(fnameR,"r");
    /*--------------------------------------------------------------------------------------*/
    /*書出用コメント入力*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1);
printf("*comment\n");
    /*--------------------------------------------------------------------------------------*/
    /*節点数，要素数，材料種別数][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,","));
    nod=Nnod;
    nhen=Nfre;
    nt=NODT*nhen;
printf("*basic data\n");
    /*------------------------------------------------------------------------------------------*/
    /*材料特性（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,","));
        wm[i][4]=atof(strtok(NULL,","));
        wm[i][5]=atof(strtok(NULL,","));
        wm[i][6]=atof(strtok(NULL,","));
        wm[i][7]=atof(strtok(NULL,","));
    }
printf("*material\n");
    /*------------------------------------------------------------------------------------------*/
    /*要素構成節点と材料種別（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];
                gamma[ne]=wm[i][4];
                gkh[ne]  =wm[i][5];
                gkv[ne]  =wm[i][6];
                alpha[ne]=wm[i][7];
            }
        }
    }
printf("*element\n");
    /*--------------------------------------------------------------------------------------*/
    /*節点座標（x,y）*/
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=NODT;i++){
        fgets(dat,sizeof dat,fin);
        x[i]     =atof(strtok(dat,","));
        y[i]     =atof(strtok(NULL,","));
        deltaT[i]=atof(strtok(NULL,","));
    }
printf("*coordinate\n");

    /*--------------------------------------------------------------------------------------*/
    /*変位既知節点番号，既知節点変位量*/
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=nt;i++){rdis[i]=0.0;}
    if(0<KOX){
        for(i=1;i<=KOX;i++){
            fgets(dat,sizeof dat,fin);
            nokx[i]=atoi(strtok(dat,","));
            rdis[3*nokx[i]-2]=atof(strtok(NULL,","));
        }
    }
    if(0<KOY){
        for(i=1;i<=KOY;i++){
            fgets(dat,sizeof dat,fin);
            noky[i]=atoi(strtok(dat,","));
            rdis[3*noky[i]-1]=atof(strtok(NULL,","));
        }
    }
    if(0<KOZ){
        for(i=1;i<=KOZ;i++){
            fgets(dat,sizeof dat,fin);
            nokz[i]=atoi(strtok(dat,","));
            rdis[3*nokz[i]]=atof(strtok(NULL,","));
        }
    }
printf("*fixed condition\n");
    /*--------------------------------------------------------------------------------------*/
    /*節点番号,x方向荷重,y方向荷重,モーメント荷重*/
    /*--------------------------------------------------------------------------------------*/
    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[3*n-2]=atof(strtok(NULL,","));
            f[3*n-1]=atof(strtok(NULL,","));
            f[3*n]  =atof(strtok(NULL,","));
        }
    }
printf("*load\n");
    fclose(fin);

    /*******************************************************/
    /*入力確認出力*/
    /*******************************************************/
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s\n",strcom);
    fprintf(fout,"NODT,NELT,MATEL,KOX,KOY,KOZ,NF\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",KOX)   ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOY)   ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOZ)   ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",NF)    ;strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
    fprintf(fout,"*node characteristics\n");
    fprintf(fout,"node,x,y,fx,fy,fz,fix-x,fix-y,fix-z,rdis-x,rdis-y,rdis-z,deltaT\n");
    for(i=1;i<=NODT;i++){
        k=0;l=0;m=0;
        if(0<KOX){
            for(j=1;j<=KOX;j++){
                if(i==nokx[j])k=1;
            }
        }
        if(0<KOY){
            for(j=1;j<=KOY;j++){
                if(i==noky[j])l=1;
            }
        }
        if(0<KOZ){
            for(j=1;j<=KOZ;j++){
                if(i==nokz[j])m=1;
            }
        }
        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",f[3*i-2])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",f[3*i-1])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",f[3*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,"%d",m)           ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",rdis[3*i-2]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",rdis[3*i-1]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",rdis[3*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,E,A,I,gamma,gkh,gkv,alpha,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,"%e",Em[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",AA[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",AI[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",gamma[ne])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",gkh[ne])      ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",gkv[ne])      ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",alpha[ne])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",matno[ne])    ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fclose(fout);
    /*******************************************************/
    /*メイン処理計算*/
    /*******************************************************/
    /*-----------------------------------------------------*/
    /*全体剛性行列・全体荷重ベクトルクリア*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        ftvec[i]=0.0;
        fmass[i]=0.0;
        ftemp[i]=0.0;
        dis[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_FRM(ne, kakom, x, y, Em, AA, AI, sm, hen);
        ASMAT(ne, nod, nhen, kakom, sm, tsm);

        /*慣性力ベクトル組立*/
        MV_FRM(ne,kakom,x,y,AA,gamma,gkh,gkv,fevec,hen);
        ASVEC(ne, nod, nhen, kakom, fevec, fmass);

        /*温度荷重相当ベクトル組立*/
        TV_FRM(ne,kakom, deltaT, Em, AA,alpha,fevec,hen);
        ASVEC(ne, nod, nhen, kakom, fevec, ftemp);
    }
    /*-----------------------------------------------------*/
    /*既知節点変位による荷重成分ベクトル作成*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        fdis[i]=0.0;
        for(j=1;j<=nt;j++){
            fdis[i]=fdis[i]+rdis[j]*tsm[i][j];
        }
    }
    /*-----------------------------------------------------*/
    /*荷重ベクトルセット*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        ftvec[i]=f[i]+fmass[i]+ftemp[i]-fdis[i];
    }
    /*-----------------------------------------------------*/
    /*変位既知節点自由度記憶*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        index[i]=i;
    }
    if(0<KOX){
        for(i=1;i<=KOX;i++){n=3*nokx[i]-2;index[n]=0;}
    }
    if(0<KOY){
        for(i=1;i<=KOY;i++){n=3*noky[i]-1;index[n]=0;}
    }
    if(0<KOZ){
        for(i=1;i<=KOZ;i++){n=3*nokz[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];
        reac[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){
            printf("対角項<=0","異常終了通知");
            return 0;
        }
    }
    /*-----------------------------------------------------*/
    /*連立一次方程式（コレスキー法）*/
    /*-----------------------------------------------------*/
    CHBAND10(mm,ib,tsm);/*三角行列*/
    CHBAND11(mm,ib,tsm,reac);/*前進消去・後退代入*/
    /*-----------------------------------------------------*/
    /*変位計算*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){dis[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;
        dis[i]=dis[i]+rdis[i];/*既知節点変位の足し込み*/
    }
    /*-----------------------------------------------------*/
    /*断面力計算*/
    /*-----------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        SM_FRM(ne,kakom,x,y,Em,AA,AI,sm,hen);
        CALST_FRM(ne,kakom,dis,sm,hen,reac,tsm,Em,AA,alpha,deltaT);
    }
    /*******************************************************/
    /*結果出力*/
    /*******************************************************/
    fout=fopen(fnameW,"a+");
    fprintf(fout,"*displacements and forces\n");
    fprintf(fout,"node,coord-x,coord-y,dir-x,dir-y,dir-z,reac-x,reac-y,reac-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",reac[3*i-2]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",reac[3*i-1]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",reac[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 resultants\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",tsm[ne][i]);strcat(dat,",");strcat(dat,sv);
        }
        fprintf(fout,"%s\n",dat);
    }
    /*計算完了*/
    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,"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 SM_FRM(int ne,int kakom[NEmax+1][Nnod+1],double x[NDmax+1],double y[NDmax+1],
    double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double sm[Nfel+1][Nfel+1],double hen[Nfel+1][Nfel+1])
{
    /*要素剛性行列・座標変換行列*/
    int i,j,k;
    double w[Nfel+1][Nfel+1];
    double AL,EA,EI,cs,sn;

    /*座標変換行列*/
    i=kakom[ne][1];
    j=kakom[ne][2];
    AL=sqrt((x[j]-x[i])*(x[j]-x[i])+(y[j]-y[i])*(y[j]-y[i]));
    cs=(x[j]-x[i])/AL;
    sn=(y[j]-y[i])/AL;
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            sm[i][j]=0.0;
            hen[i][j]=0.0;
        }
    }
    hen[1][1]=cs;
    hen[1][2]=sn;
    hen[2][1]=-sn;
    hen[2][2]=cs;
    hen[3][3]=1.0;
    for(i=1;i<=3;i++){
        for(j=1;j<=3;j++){
            hen[i+3][j+3]=hen[i][j];
        }
    }

    /*要素剛性行列*/
    EA=Em[ne]*AA[ne];
    EI=Em[ne]*AI[ne];
    sm[1][1]=EA/AL;
    sm[2][2]=12.0*EI/(AL*AL*AL);
    sm[3][2]=6.0*EI/(AL*AL);
    sm[3][3]=4.0*EI/AL;
    sm[4][1]=-EA/AL;
    sm[4][4]=sm[1][1];
    sm[5][2]=-sm[2][2];
    sm[5][3]=-sm[3][2];
    sm[5][5]=sm[2][2];
    sm[6][2]=sm[3][2];
    sm[6][3]=2.0*EI/AL;
    sm[6][5]=-sm[3][2];
    sm[6][6]=sm[3][3];
    /*対称項セット*/
    for(i=1;i<=6;i++){
        for(j=1;j<=i;j++){
            sm[j][i]=sm[i][j];
        }
    }
    /*[w]=[T]/**[k]*/
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            w[i][j]=0.0;
            for(k=1;k<=6;k++){
                w[i][j]=w[i][j]+hen[k][i]*sm[k][j];
            }
        }
    }
    /*[w]*[T]=[T]/*[k][T]：[sm]は全体座標系での要素剛性行列*/
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            sm[i][j]=0.0;
            for(k=1;k<=6;k++){
                sm[i][j]=sm[i][j]+w[i][k]*hen[k][j];
            }
        }
    }
}
/*--------------------------------------------------------------------------------*/

void MV_FRM(int ne,int kakom[NEmax+1][Nnod+1],double x[NDmax+1],double y[NDmax+1],
    double AA[NEmax+1],double gamma[NEmax+1],double gkh[NEmax+1],double gkv[NEmax+1],
    double fevec[Nfel+1],double hen[Nfel+1][Nfel+1])
{
    int i,j,k;
    double AL,am[7][7],w[7][7],work[7];
    double ggg;

    ggg=9.8;

    /*行列初期化*/
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            am[i][j]=0.0;
        }
    }

    /*要素長さ*/
    i=kakom[ne][1];
    j=kakom[ne][2];
    AL=sqrt((x[j]-x[i])*(x[j]-x[i])+(y[j]-y[i])*(y[j]-y[i]));

    /*質量行列*/
    am[1][1]=1.0/3.0;
    am[2][1]=0.0;
    am[3][1]=0.0;
    am[4][1]=1.0/6.0;
    am[5][1]=0.0;
    am[6][1]=0.0;
    am[2][2]=13.0/35.0;
    am[3][2]=11.0*AL/210.0;
    am[4][2]=0.0;
    am[5][2]=9.0/70.0;
    am[6][2]=-13.0*AL/420.0;
    am[3][3]=AL*AL/105.0;
    am[4][3]=0.0;
    am[5][3]=13.0*AL/420.0;
    am[6][3]=-AL*AL/140.0;
    am[4][4]=1.0/3.0;
    am[5][4]=0.0;
    am[6][4]=0.0;
    am[5][5]=13.0/35.0;
    am[6][5]=-11.0*AL/210.0;
    am[6][6]=AL*AL/105.0;
    /*対称項セット*/
    for(i=1;i<=5;i++){
        for(j=i+1;j<=6;j++){
            am[i][j]=am[j][i];
        }
    }
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            am[i][j]=gamma[ne]*AA[ne]*AL*am[i][j]/ggg;
        }
    }

    /*全体座標系への変換*/
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            w[i][j]=0.0;
            for(k=1;k<=6;k++){
                w[i][j]=w[i][j]+hen[k][i]*am[k][j];
            }
        }
    }
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            am[i][j]=0.0;
            for(k=1;k<=6;k++){
                am[i][j]=am[i][j]+w[i][k]*hen[k][j];
            }
        }
    }

    work[1]=ggg*gkh[ne];
    work[2]=ggg*gkv[ne];
    work[3]=0.0;
    work[4]=ggg*gkh[ne];
    work[5]=ggg*gkv[ne];
    work[6]=0.0;

    /*全体座標系での慣性力ベクトル*/
    for(i=1;i<=6;i++){
        fevec[i]=0.0;
        for(j=1;j<=6;j++){
            fevec[i]=fevec[i]+am[i][j]*work[j];
        }
    }
}
/*--------------------------------------------------------------------------------*/

void TV_FRM(int ne,int kakom[NEmax+1][Nnod+1],double deltaT[NDmax+1],
    double Em[NEmax+1],double AA[NEmax+1],double alpha[NEmax+1],double fevec[Nfel+1],double hen[Nfel+1][Nfel+1])
{
    /*温度荷重相当ベクトル*/
    int i, j;
    double tempa;
    double work[Nfel+1];

    i = kakom[ne][1];
    j = kakom[ne][2];
    tempa=0.5*(deltaT[i]+deltaT[j]);
    /*要素座標系での温度荷重相当ベクトル*/
    work[1] = -1.0*Em[ne] * AA[ne] * alpha[ne] * tempa;
    work[2] = 0.0;
    work[3] = 0.0;
    work[4] = 1.0 * Em[ne] * AA[ne] * alpha[ne] * tempa;
    work[5] = 0.0;
    work[6] = 0.0;
    /*全体座標系での温度荷重相当ベクトル*/
    /*（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる）*/
    for (i = 1; i <= 6; i++){
        fevec[i] = 0.0;
        for (j = 1; j <= 6; j++){
            fevec[i] = fevec[i] + hen[j][i] * work[j];
        }
    }
}
/*--------------------------------------------------------------------------------*/

void CALST_FRM(int ne, int kakom[NEmax+1][Nnod+1], double dis[NTmax+1], double sm[Nfel+1][Nfel+1], double hen[Nfel+1][Nfel+1], double reac[NTmax+1],
    double tsm[NTmax+1][NTmax+1],double Em[NEmax+1],double AA[NEmax+1],double alpha[NEmax+1],double deltaT[NDmax+1])
{
    /*断面力算定*/
    int i,j,k;
    double wd[Nfel+1];
    double ws[Nfel+1];
    int ia,ir,iw;
    double tempa;

    /*要素節点変位取得（全体座標系）*/
    for( i = 1 ;i<= 2;i++){
        ia = kakom[ne][i];
        wd[3 * i - 2] = dis[3 * ia - 2];
        wd[3 * i - 1] = dis[3 * ia - 1];
        wd[3 * i] = dis[3 * ia];
    }
    /*全体座標系での要素節点力算定*/
    for( i = 1 ;i<= 6;i++){
        ws[i] = 0.0;
        for( k = 1 ;k<= 6;k++){
            ws[i] = ws[i] + sm[i][k] * wd[k];
    }}
    /*要素断面力算定：座標変換行列×全体座標系要素節点力*/
    /*今後使用しない全体剛性行列格納配列[tsm]に断面力を格納*/
    for( i = 1 ;i<= 6;i++){
        tsm[ne][i] = 0.0;
        for( k = 1 ;k<= 6;k++){
            tsm[ne][i] = tsm[ne][i] + hen[i][k] * ws[k];
    }}
    /*軸力の温度応力補正*/
    i = kakom[ne][1];
    j = kakom[ne][2];
    tempa = 0.5 * (deltaT[i] + deltaT[j]);
    tsm[ne][1] = tsm[ne][1] + Em[ne] * AA[ne] * alpha[ne] * tempa;
    tsm[ne][4] = tsm[ne][4] - Em[ne] * AA[ne] * alpha[ne] * tempa;
    /*全体座標系での節点荷重ベクトル*/
    /*（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる）*/
    for (i = 1; i <= 6; i++){
        ws[i] = 0.0;
        for (j = 1; j <= 6; j++){
            ws[i] = ws[i] + hen[j][i] * tsm[ne][j];
        }
    }
    /*要素節点力の全体系へのセット*/
    for (i = 1; i <= 2; i++){
        ia = kakom[ne][i];
        for (j = 1; j <= 3; j++){
            ir = 3 * (ia - 1) + j;
            iw = 3 * (i - 1) + j;
            reac[ir] = reac[ir] + ws[iw];
        }
    }
}
/*--------------------------------------------------------------------------------*/

void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][Nnod+1],double sm[Nfel+1][Nfel+1],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][Nnod+1],double fevec[Nfel+1],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];
        }
    }
}
/*--------------------------------------------------------------------------------*/