/*===========================*/
/* gccMESH4nod-pg.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

int main_part(char fnameR[50],char fnameF[50]);
void GPL_mesh(char fnameF[50],int NODT,int NELT,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1]);
void GPL_stress(char fnameF[50],int NODT,int NELT,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],double stress[NEmax+1][5][7]);
/*------------------------------------------------------------------------------*/

int main(int argc, char* argv[])
{
    char fnameR[50];
    char fnameF[50];

    strcpy(fnameR,argv[1]);
    strcpy(fnameF,argv[2]);

    main_part(fnameR,fnameF);

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

int main_part(char fnameR[50],char fnameF[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 NSTRES;           /* 応力状態（0;平面歪み,1;平面応力） */
    int NODT;             /* 節点総数 */
    int NELT;             /* 要素総数 */
    int MATEL;            /* 材料種類数 */
    int KOX;              /* x方向変位既知節点数 */
    int KOY;              /* y方向変位既知節点数 */
    int NF;               /* 載荷節点数 */
    /*--------------------------------------------------------------------------------------*/
    /*配列寸法宣言*/
    /*--------------------------------------------------------------------------------------*/
    int    kakom[NEmax+1][5];         /* 要素構成節点番号 */
    int    matno[NEmax+1];            /* 材料種別No */
    double x[NDmax+1];                /* 節点x座標 */
    double y[NDmax+1];                /* 節点y座標 */
    double wk;
    int iw,iw1,iw2;
    double stress[NEmax+1][5][7];

    /******************************************************/
    /*データ入力*/
    /******************************************************/
    fin=fopen(fnameR,"r");
    fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1);
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);
    NSTRES=atoi(strtok(dat,","));
    NODT  =atoi(strtok(NULL,","));
    NELT  =atoi(strtok(NULL,","));
    MATEL =atoi(strtok(NULL,","));
    KOX   =atoi(strtok(NULL,","));
    KOY   =atoi(strtok(NULL,","));
    NF    =atoi(strtok(NULL,","));
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);
    for(i=1;i<=NODT;i++){
        fgets(dat,sizeof dat,fin);
        iw  =atoi(strtok(dat,","));
        x[i]=atof(strtok(NULL,","));
        y[i]=atof(strtok(NULL,","));
    }
    fgets(dat,sizeof dat,fin);
    for(ne=1;ne<=NELT;ne++){
        fgets(dat,sizeof dat,fin);
        iw=atoi(strtok(dat,","));
        kakom[ne][1]=atoi(strtok(NULL,","));
        kakom[ne][2]=atoi(strtok(NULL,","));
        kakom[ne][3]=atoi(strtok(NULL,","));
        kakom[ne][4]=atoi(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        wk=atof(strtok(NULL,","));
        matno[ne]=atoi(strtok(NULL,","));
    }
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);
    for(i=1;i<=NODT;i++){
        fgets(dat,sizeof dat,fin);
    }
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);
    for(i=1;i<=NODT;i++){
        fgets(dat,sizeof dat,fin);
    }
    fgets(dat,sizeof dat,fin);
    fgets(dat,sizeof dat,fin);
    for(ne=1;ne<=NELT;ne++){
        fgets(dat,sizeof dat,fin);
            iw1=atoi(strtok(dat,","));
            iw2=atoi(strtok(NULL,","));
            stress[iw1][iw2][1]=atof(strtok(NULL,","));
            stress[iw1][iw2][2]=atof(strtok(NULL,","));
            stress[iw1][iw2][3]=atof(strtok(NULL,","));
            stress[iw1][iw2][4]=atof(strtok(NULL,","));
            stress[iw1][iw2][5]=atof(strtok(NULL,","));
            stress[iw1][iw2][6]=atof(strtok(NULL,","));
        fgets(dat,sizeof dat,fin);
            iw1=atoi(strtok(dat,","));
            iw2=atoi(strtok(NULL,","));
            stress[iw1][iw2][1]=atof(strtok(NULL,","));
            stress[iw1][iw2][2]=atof(strtok(NULL,","));
            stress[iw1][iw2][3]=atof(strtok(NULL,","));
            stress[iw1][iw2][4]=atof(strtok(NULL,","));
            stress[iw1][iw2][5]=atof(strtok(NULL,","));
            stress[iw1][iw2][6]=atof(strtok(NULL,","));
        fgets(dat,sizeof dat,fin);
            iw1=atoi(strtok(dat,","));
            iw2=atoi(strtok(NULL,","));
            stress[iw1][iw2][1]=atof(strtok(NULL,","));
            stress[iw1][iw2][2]=atof(strtok(NULL,","));
            stress[iw1][iw2][3]=atof(strtok(NULL,","));
            stress[iw1][iw2][4]=atof(strtok(NULL,","));
            stress[iw1][iw2][5]=atof(strtok(NULL,","));
            stress[iw1][iw2][6]=atof(strtok(NULL,","));
        fgets(dat,sizeof dat,fin);
            iw1=atoi(strtok(dat,","));
            iw2=atoi(strtok(NULL,","));
            stress[iw1][iw2][1]=atof(strtok(NULL,","));
            stress[iw1][iw2][2]=atof(strtok(NULL,","));
            stress[iw1][iw2][3]=atof(strtok(NULL,","));
            stress[iw1][iw2][4]=atof(strtok(NULL,","));
            stress[iw1][iw2][5]=atof(strtok(NULL,","));
            stress[iw1][iw2][6]=atof(strtok(NULL,","));
    }
    fclose(fin);
/*
    GPL_mesh(fnameF,NODT,NELT,kakom,x,y);
*/
    GPL_stress(fnameF,NODT,NELT,kakom,x,y,stress);
}
/*------------------------------------------------------------------------------*/
void GPL_mesh(char fnameF[50],int NODT,int NELT,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1])
{
    FILE *fp;
    char fnameW[50]="gpl_temp.txt";
    double x0,y0, x1,y1,x2,y2,x3,y3,ang;
    char str[50],sv[50],cn[10];
    int i,j,k,l,ne;
    double xmin,xmax,ymin,ymax;

    strcpy(cn,"#FFFFE0");
    xmin=x[1];
    xmax=x[1];
    ymin=y[1];
    ymax=y[1];
    for(i=1;i<=NODT;i++){
        if(x[i]<xmin)xmin=x[i];
        if(xmax<x[i])xmax=x[i];
        if(y[i]<ymin)ymin=y[i];
        if(ymax<y[i])ymax=y[i];
    }

    fp=fopen(fnameW,"w");
    fprintf(fp,"reset\n");
    fprintf(fp,"set terminal postscript eps enhanced color font \"Arial\" 12\n");
    fprintf(fp,"set output \"%s\"\n",fnameF);
    fprintf(fp,"#\n");
    fprintf(fp,"set size ratio -1\n");
    fprintf(fp,"#\n");
    fprintf(fp,"set xrange [%.3f:%.3f]\n",xmin,xmax);
    fprintf(fp,"set yrange [%.3f:%.3f]\n",ymin,ymax);
    fprintf(fp,"set border 0\n");
    fprintf(fp,"set xtics 0,0,0\n");
    fprintf(fp,"set ytics 0,0,0\n");
    fprintf(fp,"set xtics nomirror\n");
    fprintf(fp,"set ytics nomirror\n");
    fprintf(fp,"#\n");

    fprintf(fp,"#\n");

    for(ne=1;ne<=NELT;ne++){
        i=kakom[ne][1];
        j=kakom[ne][2];
        k=kakom[ne][3];
        l=kakom[ne][4];
        fprintf(fp,"set object %d polygon from %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f\n",
            ne,x[i],y[i],x[j],y[j],x[k],y[k],x[l],y[l],x[i],y[i]);
        fprintf(fp,"set object %d fillcolor rgb \"%s\" fillstyle solid 1.0 border lt -1\n",ne,cn);
    }
    fprintf(fp,"#\n");
    fprintf(fp,"plot \"-\" with point notitle\n");
    fprintf(fp,"%.3f %.3f\n",x[1],y[1]);
    fprintf(fp,"e\n");

    fprintf(fp,"#\n");
    fprintf(fp,"quit\n");
    fclose(fp);

    fp=popen("gnuplot.exe","w");
        fprintf(fp,"load \"%s\"\n",fnameW);
        fflush(fp);
    pclose(fp);

}
/*---------------------------------------------------------------------------*/
void GPL_stress(char fnameF[50],int NODT,int NELT,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],double stress[NEmax+1][5][7])
{
    FILE *fp;
    char fnameW[50]="gpl_temp.txt";
    char str[50],sv[50],cn[10];
    int i,j,k,l,ne,kk,jj;
    double xmin,xmax,ymin,ymax;
    double pmax,pmin,ss;
    double x11,x22,x33,x44,x12,x23,x34,x41,x00;
    double y11,y22,y33,y44,y12,y23,y34,y41,y00;

    pmax=stress[1][1][4];
    pmin=stress[1][1][4];
    for(ne=1;ne<=NELT;ne++){
        for(kk=1;kk<=4;kk++){
            if(stress[ne][kk][4]<pmin)pmin=stress[ne][kk][4];
            if(pmax<stress[ne][kk][4])pmax=stress[ne][kk][4];
            if(stress[ne][kk][5]<pmin)pmin=stress[ne][kk][5];
            if(pmax<stress[ne][kk][5])pmax=stress[ne][kk][5];
    }}

    xmin=x[1];
    xmax=x[1];
    ymin=y[1];
    ymax=y[1];
    for(i=1;i<=NODT;i++){
        if(x[i]<xmin)xmin=x[i];
        if(xmax<x[i])xmax=x[i];
        if(y[i]<ymin)ymin=y[i];
        if(ymax<y[i])ymax=y[i];
    }

    fp=fopen(fnameW,"w");
    fprintf(fp,"reset\n");
    fprintf(fp,"set terminal postscript eps enhanced color font \"Arial\" 12\n");
    fprintf(fp,"set output \"%s\"\n",fnameF);
    fprintf(fp,"#\n");
    fprintf(fp,"set size ratio -1\n");
    fprintf(fp,"#\n");
    fprintf(fp,"set xrange [%.3f:%.3f]\n",xmin,xmax);
    fprintf(fp,"set yrange [%.3f:%.3f]\n",ymin,ymax);
    fprintf(fp,"set border 0\n");
    fprintf(fp,"set xtics 0,0,0\n");
    fprintf(fp,"set ytics 0,0,0\n");
    fprintf(fp,"set xtics nomirror\n");
    fprintf(fp,"set ytics nomirror\n");
    fprintf(fp,"#\n");
    fprintf(fp,"set rmargin 11\n");

    fprintf(fp,"#\n");
    fprintf(fp,"set multiplot layout 1,2\n");

    /*jj=4;第一主応力*/
    /*jj=5;第二主応力*/
    for(jj=4;jj<=5;jj++){
        for(ne=1;ne<=NELT;ne++){
            i=kakom[ne][1];
            j=kakom[ne][2];
            k=kakom[ne][3];
            l=kakom[ne][4];
            x11=x[i];y11=y[i];
            x22=x[j];y22=y[j];
            x33=x[k];y33=y[k];
            x44=x[l];y44=y[l];
            x12=0.5*(x11+x22);y12=0.5*(y11+y22);
            x23=0.5*(x22+x33);y23=0.5*(y22+y33);
            x34=0.5*(x33+x44);y34=0.5*(y33+y44);
            x41=0.5*(x44+x11);y41=0.5*(y44+y11);
            x00=0.5*(x11+x33);y00=0.5*(y11+y33);
            for(kk=1;kk<=4;kk++){
                ss=stress[ne][kk][jj];
                if(0.0<=pmax||0.0<pmin){
                    if(pmax*5.0/6.0<=ss)strcpy(cn,"#0000ff");
                    if(pmax*4.0/6.0<=ss&&ss<pmax*5.0/6.0)strcpy(cn,"#0033ff");
                    if(pmax*3.0/6.0<=ss&&ss<pmax*4.0/6.0)strcpy(cn,"#0066ff");
                    if(pmax*2.0/6.0<=ss&&ss<pmax*3.0/6.0)strcpy(cn,"#0099ff");
                    if(pmax*1.0/6.0<=ss&&ss<pmax*2.0/6.0)strcpy(cn,"#00ccff");
                    if(pmax*0.0/6.0<=ss&&ss<pmax*1.0/6.0)strcpy(cn,"#00ffff");
                }
                if(pmin<=0.0||pmax<0.0){
                    if(pmin*1.0/6.0<=ss&&ss<pmin*0.0/6.0)strcpy(cn,"#ffff00");
                    if(pmin*2.0/6.0<=ss&&ss<pmin*1.0/6.0)strcpy(cn,"#ffcc00");
                    if(pmin*3.0/6.0<=ss&&ss<pmin*2.0/6.0)strcpy(cn,"#ff9900");
                    if(pmin*4.0/6.0<=ss&&ss<pmin*3.0/6.0)strcpy(cn,"#ff6600");
                    if(pmin*5.0/6.0<=ss&&ss<pmin*4.0/6.0)strcpy(cn,"#ff3300");
                    if(ss<pmin*5.0/6.0)strcpy(cn,"#ff0000");
                }
                switch(kk){
                    case 1:
                        fprintf(fp,"set object %d polygon from %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f\n",
                            4*ne-3,x11,y11,x12,y12,x00,y00,x41,y41,x11,y11);
                    break;
                    case 2:
                        fprintf(fp,"set object %d polygon from %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f\n",
                            4*ne-2,x22,y22,x23,y23,x00,y00,x12,y12,x22,y22);
                    break;
                    case 3:
                        fprintf(fp,"set object %d polygon from %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f\n",
                            4*ne-1,x33,y33,x34,y34,x00,y00,x23,y23,x33,y33);
                    break;
                    case 4:
                        fprintf(fp,"set object %d polygon from %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f to %.3f,%.3f\n",
                            4*ne-0,x44,y44,x41,y41,x00,y00,x34,y34,x44,y44);
                    break;
                }
                fprintf(fp,"set object %d fillcolor rgb \"%s\" fillstyle solid 1.0\n",4*(ne-1)+kk,cn);
            }
        }
        fprintf(fp,"#\n");
        i=4*NELT;
        if(0.0<=pmax||0.0<pmin){
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.95 to graph 1.1,1.00 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#0000ff");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.90 to graph 1.1,0.95 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#0033ff");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.85 to graph 1.1,0.90 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#0066ff");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.80 to graph 1.1,0.85 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#0099ff");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.75 to graph 1.1,0.80 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#00ccff");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.70 to graph 1.1,0.75 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#00ffff");
            fprintf(fp,"set label 1 \" %.3f\" at graph 1.10,graph 1.00 left\n",pmax);
            fprintf(fp,"set label 2 \" %.3f\" at graph 1.10,graph 0.70 left\n",0.0);
        }
        if(pmin<=0.0||pmax<0.0){
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.65 to graph 1.1,0.70 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#ffff00");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.60 to graph 1.1,0.65 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#ffcc00");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.55 to graph 1.1,0.60 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#ff9900");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.50 to graph 1.1,0.55 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#ff6600");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.45 to graph 1.1,0.50 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#ff3300");
            i=i+1;fprintf(fp,"set object %d rect from graph 1.05,0.40 to graph 1.1,0.45 fillcolor rgb \"%s\" fillstyle solid 1.0 border -1\n",i,"#ff0000");
            fprintf(fp,"set label 3 \" %.3f\" at graph 1.10,graph 0.70 left\n",0.0);
            fprintf(fp,"set label 4 \" %.3f\" at graph 1.10,graph 0.40 left\n",pmin);
        }
        fprintf(fp,"#\n");
        fprintf(fp,"plot \"-\" with dots notitle\n");
        fprintf(fp,"%.3f %.3f\n",xmin,ymin);
        fprintf(fp,"%.3f %.3f\n",xmax,ymin);
        fprintf(fp,"%.3f %.3f\n",xmax,ymax);
        fprintf(fp,"%.3f %.3f\n",xmin,ymax);
        fprintf(fp,"%.3f %.3f\n",xmin,ymin);
        fprintf(fp,"e\n");
    }

    fprintf(fp,"unset multiplot\n");
    fprintf(fp,"#\n");
    fprintf(fp,"quit\n");
    fclose(fp);

    fp=popen("gnuplot.exe","w");
        fprintf(fp,"load \"%s\"\n",fnameW);
        fflush(fp);
    pclose(fp);

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