#include<stdio.h>
#include<math.h>
main()
{
   int np,np2,ne,p,koteix[100],koteiy[100],i,j,k,l,ielm[100][3];
   float x[100],y[100],em[6][6],tk[150][150],loadx[100],loady[100],f[150];
   float d[3][3],bb[3][6],bbt[6][3],btd[6][3];
   float eps[100][3],sigma[100][3],disp[150];
   float elas,pois,s,b[3],c[3],c1,c2,pmax,w,ww[150];
   float x0,x1,x2,y0,y1,y2,z;
   FILE *fp1;
   FILE *fp2;
/* ----- zero clear of matrix and vector ------------- */
   for(i=0;i<100;i++){
     x[i]=0.0;
     y[i]=0.0;
     koteix[i]=0;
	 koteiy[i]=0;
     loadx[i]=0.0;
     loady[i]=0.0;
   }
   for(i=0;i<100;i++){
     ielm[i][0]=0;
     ielm[i][1]=0;
     ielm[i][2]=0;
     eps[i][0]=0.0;
     eps[i][1]=0.0;
     eps[i][2]=0.0;
     sigma[i][0]=0.0;
     sigma[i][1]=0.0;
     sigma[i][2]=0.0;
   }
   for(i=0;i<150;i++){
     disp[i]=0.0;
     f[i]=0.0;
     ww[i]=0.0;
     for(j=0;j<150;j++){
       tk[i][j]=0.0;
     }
   }
   for(i=0;i<3;i++){
     for(j=0;j<3;j++){
       d[i][j]=0.0;
     }
   }
/* ------------ input of data ------------------------ */
   if((fp1=fopen("c:data3.txt","r"))==NULL){
     printf("no data file\n");
	 exit(1);
   }
     fscanf(fp1,"%d %d",&np,&ne);
     fscanf(fp1,"%f %f",&elas,&pois);
     for(i=0;i<np;i++){
       fscanf(fp1,"%f %f",&x[i],&y[i]);
     }
     for(i=0;i<ne;i++){
       fscanf(fp1,"%d %d %d",&ielm[i][0],&ielm[i][1],&ielm[i][2]);
     }
     for(i=0;i<np;i++){
       fscanf(fp1,"%d %d",&koteix[i],&koteiy[i]);
     }
     for(i=0;i<np;i++){
       fscanf(fp1,"%f %f",&loadx[i],&loady[i]);
     }   
     fclose(fp1);
/* ----- input data check ----- */	
	 printf("%5d %5d\n",np,ne);
     printf("%f %f\n",elas,pois);
     for(i=0;i<np;i++){
       printf("%f %f\n",x[i],y[i]);
     }
     for(i=0;i<ne;i++){
       printf("%5d %5d %5d\n",ielm[i][0],ielm[i][1],ielm[i][2]);
     }
     for(i=0;i<np;i++){
       printf("%d %d\n",koteix[i],koteiy[i]);
     }
     for(i=0;i<np;i++){
       printf("%f %f\n",loadx[i],loady[i]);
     }  
/* ------  d matrix -------- */
   c1=elas/(1.0-pois*pois);
   c2=pois;
   d[0][0]=c1;
   d[0][1]=c1*c2;
   d[1][0]=c1*c2;
   d[1][1]=c1;
   d[2][2]=elas/(2.0*(1.0+pois));
/* ----- output of d-matrix ----- */
   for(i=0;i<3;i++){
	   for(j=0;j<3;j++){
		   printf("%10.2e",d[i][j]);
       }
	   printf("\n");
   }
/* ----------------------------- */
   np2=np*2;
/* -----  element matrix ----- */
   for(i=0;i<ne;i++){
/* ---------- zero clear of matrix ------------ */
     for(j=0;j<6;j++){
       for(k=0;k<6;k++){
         em[j][k]=0.0;
       }
     }
     for(j=0;j<3;j++){
       b[j]=0.0;
       c[j]=0.0;
     }
     for(j=0;j<6;j++){
       for(k=0;k<3;k++){
         bb[k][j]=0.0;
         bbt[j][k]=0.0;
         btd[j][k]=0.0;
     }
   }
/* ---------------------------------------------- */
   x0=x[ielm[i][0]-1];
   x1=x[ielm[i][1]-1];
   x2=x[ielm[i][2]-1];
   y0=y[ielm[i][0]-1];
   y1=y[ielm[i][1]-1];
   y2=y[ielm[i][2]-1];
/* -------------------- */
   b[0]=y1-y2;
   b[1]=y2-y0;
   b[2]=y0-y1;
   c[0]=x2-x1;
   c[1]=x0-x2;
   c[2]=x1-x0;
   s=fabs(x1*y2+x2*y0+x0*y1-x1*y0-x0*y2-x2*y1)/2.0;
/* -------------------------- */   
   bb[0][0]=b[0]/(2.0*s);
   bb[0][2]=b[1]/(2.0*s);
   bb[0][4]=b[2]/(2.0*s);
   bb[1][1]=c[0]/(2.0*s);
   bb[1][3]=c[1]/(2.0*s);
   bb[1][5]=c[2]/(2.0*s);
   bb[2][0]=c[0]/(2.0*s);
   bb[2][1]=b[0]/(2.0*s);
   bb[2][2]=c[1]/(2.0*s);
   bb[2][3]=b[1]/(2.0*s);
   bb[2][4]=c[2]/(2.0*s);
   bb[2][5]=b[2]/(2.0*s);
/* ------------------------ */
   for(j=0;j<6;j++){
     for(k=0;k<3;k++){
       bbt[j][k]=bb[k][j];
     }
   }
/* ----- output of b-matrix and bt-matrix ----- */
   printf("---------------------\n");
   printf("i=%d \n",i+1);
   printf("x0=%f y0=%f\n",x0,y0);
   printf("x1=%f y1=%f\n",x1,y1);
   printf("x2=%f y2=%f\n",x2,y2);
   printf("s=%f\n",s);
   for(j=0;j<3;j++){
     for(k=0;k<6;k++){
       printf("%10.2e",bb[j][k]);
     }
	 printf("\n");
   }
   printf("\n");
   for(j=0;j<6;j++){
     for(k=0;k<3;k++){
       printf("%10.2e",bbt[j][k]);
     }
	 printf("\n");
   }
   printf("---------------------\n");
/* ----------------------------- */   
   for(j=0;j<6;j++){
     for(k=0;k<3;k++){
       btd[j][k]=0.0;
       for(l=0;l<3;l++){
         btd[j][k]=btd[j][k]+bbt[j][l]*d[l][k];
       }
     }
   }
/* ----- output of btd-matrix ----- */
    for(j=0;j<6;j++){
      for(k=0;k<3;k++){
        printf("%10.2e",btd[j][k]);
      }
	  printf("\n");
   }
/* ------------------------------- */
   for(j=0;j<6;j++){
     for(k=0;k<6;k++){
       em[j][k]=0.0;
       for(l=0;l<3;l++){
         em[j][k]=em[j][k]+btd[j][l]*bb[l][k];
       }
     }
   }
/* ------------------------------------ */
  for(j=0;j<6;j++){
     for(k=0;k<6;k++){
       em[j][k]=s*1.0*em[j][k];
     }
   } 
/* ----- output of em-matrix ----- */
   printf("\n");
   for(j=0;j<6;j++){
     for(k=0;k<6;k++){
       printf("%10.2e",em[j][k]);
     }
	 printf("\n");
   }
   printf("\n");
/* ------- superposition -------------- */
   for(j=0;j<3;j++){
     for(k=0;k<3;k++){
       tk[ielm[i][j]*2-2][ielm[i][k]*2-2]=tk[ielm[i][j]*2-2][ielm[i][k]*2-2]
                                         +em[(j+1)*2-2][(k+1)*2-2];
       tk[ielm[i][j]*2-2][ielm[i][k]*2-1]=tk[ielm[i][j]*2-2][ielm[i][k]*2-1]
                                         +em[(j+1)*2-2][(k+1)*2-1];
       tk[ielm[i][j]*2-1][ielm[i][k]*2-2]=tk[ielm[i][j]*2-1][ielm[i][k]*2-2]
                                         +em[(j+1)*2-1][(k+1)*2-2];
       tk[ielm[i][j]*2-1][ielm[i][k]*2-1]=tk[ielm[i][j]*2-1][ielm[i][k]*2-1]
                                         +em[(j+1)*2-1][(k+1)*2-1];
     }
   }
 }  
/* ------ boundary condition -------- */
  for(i=0;i<np;i++){
     if(koteix[i]==1){
       for(j=0;j<np2;j++){
         tk[j][(i+1)*2-2]=0.0;
         tk[(i+1)*2-2][j]=0.0;
       }  
         tk[(i+1)*2-2][(i+1)*2-2]=1.0;
     }
  }
  for(i=0;i<np;i++){
	  if(koteiy[i]==1){
        for(j=0;j<np2;j++){
          tk[j][(i+1)*2-1]=0.0;
          tk[(i+1)*2-1][j]=0.0;
        }  
          tk[(i+1)*2-1][(i+1)*2-1]=1.0;
	  }
  }	  
/* ----------- load vector --------------- */
   for(i=0;i<np;i++){
    f[(i+1)*2-2]=loadx[i];
    f[(i+1)*2-1]=loady[i];
   }
/* ------ output tk-matrix and f-vector ----- */
   for(i=0;i<np2;i++){
	   for(j=0;j<np2;j++){
         printf("%14.2e",tk[i][j]);
       }
	   printf("\n");
   }
   printf("\n");
   for(i=0;i<np2;i++){
     printf("%10.2e",f[i]);
   }
   printf("\n");
/* ------ ガウスの消去法（ピボット選択）----- */
/* 前進消去（ピボット選択）*/

   for(k = 0; k < np2-1; k++){  /* 第ｋステップ */
      p = k;
      pmax = fabs( tk[k][k] );
      for(i = k+1; i < np2; i++){  /* ピボット選択 */
         if(fabs( tk[i][k] ) > pmax){
            p = i;
            pmax = fabs( tk[i][k] );
         }
      }

 /* エラー処理：ピボットがあまりに小さい時はメッセージを表示して終了　*/
      if(fabs( pmax ) < 1.0e-12){
         printf("too small pivot! \n");
         return(0);
      }
      if(p != k){  /* 第ｋ行と第ｐ行の交換　*/
         for(i = k; i < np2; i++){
            /* 係数行列　*/
            w = tk[k][i];
            tk[k][i] = tk[p][i];
            tk[p][i] = w;
         }
         /* 既知ベクトル */
         w = f[k];
         f[k] = f[p];
         f[p] = w;
      }
/* 前進消去 */
      for(i = k +1; i < np2; i++){ /* 第ｉ行 */
         ww[i] = tk[i][k] / tk[k][k];
         tk[i][k] = 0.0;
         /* 第ｋ行を-tk[i][k]/tk[k][k]倍して、第ｉ行に加える */
         for(j = k + 1; j < np2; j++){
            tk[i][j] = tk[i][j] - tk[k][j] * ww[i];
         }
         f[i] = f[i] - f[k] * ww[i];
      }
   }
/* 後退代入 */
      for(i = np2 - 1; i >= 0; i--){
         for(j = i + 1; j < np2; j++){
            f[i] = f[i] - tk[i][j] * f[j];
            tk[i][j] = 0.0;
         }
         f[i] = f[i] / tk[i][i];
         tk[i][i] = 1.0;
      }
      for(i=0;i<np2;i++){
        disp[i]=f[i];
      }
/* ----- output of displacement ----- */
	  for(i=0;i<np2;i++){
        printf("i=%4d %10.3e\n",i+1,disp[i]);
      }
/* -------- strain and stress of element ------------------- */
   for(i=0;i<ne;i++){
     
     for(j=0;j<3;j++){
       b[j]=0.0;
       c[j]=0.0;
     }
     for(j=0;j<6;j++){
       for(k=0;k<3;k++){
         bb[k][j]=0.0;
     }
   }
/* ---------------------------------------------- */
   x0=x[ielm[i][0]-1];
   x1=x[ielm[i][1]-1];
   x2=x[ielm[i][2]-1];
   y0=y[ielm[i][0]-1];
   y1=y[ielm[i][1]-1];
   y2=y[ielm[i][2]-1];
/* -------------------- */
   b[0]=y1-y2;
   b[1]=y2-y0;
   b[2]=y0-y1;
   c[0]=x2-x1;
   c[1]=x0-x2;
   c[2]=x1-x0;
   s=fabs(x1*y2+x2*y0+x0*y1-x1*y0-x0*y2-x2*y1)/2.0;
/* -------------------------- */   
   bb[0][0]=b[0]/(2.0*s);
   bb[0][2]=b[1]/(2.0*s);
   bb[0][4]=b[2]/(2.0*s);
   bb[1][1]=c[0]/(2.0*s);
   bb[1][3]=c[1]/(2.0*s);
   bb[1][5]=c[2]/(2.0*s);
   bb[2][0]=c[0]/(2.0*s);
   bb[2][1]=b[0]/(2.0*s);
   bb[2][2]=c[1]/(2.0*s);
   bb[2][3]=b[1]/(2.0*s);
   bb[2][4]=c[2]/(2.0*s);
   bb[2][5]=b[2]/(2.0*s);
/* --------------------------------------------------------------------- */
   eps[i][0]=bb[0][0]*disp[ielm[i][0]*2-2]+bb[0][2]*disp[ielm[i][1]*2-2]
            +bb[0][4]*disp[ielm[i][2]*2-2];
   eps[i][1]=bb[1][1]*disp[ielm[i][0]*2-1]+bb[1][3]*disp[ielm[i][1]*2-1]
            +bb[1][5]*disp[ielm[i][2]*2-1]; 
   eps[i][2]=bb[2][0]*disp[ielm[i][0]*2-2]+bb[2][1]*disp[ielm[i][0]*2-1]
            +bb[2][2]*disp[ielm[i][1]*2-2]+bb[2][3]*disp[ielm[i][1]*2-1]
            +bb[2][4]*disp[ielm[i][2]*2-2]+bb[2][5]*disp[ielm[i][2]*2-1];
/* ---------------------------------------------------------------------- */
   sigma[i][0]=d[0][0]*eps[i][0]+d[0][1]*eps[i][1];
   sigma[i][1]=d[1][0]*eps[i][0]+d[1][1]*eps[i][1];
   sigma[i][2]=d[2][2]*eps[i][2];
/* ---------------------------------------------------------------------- */
}
/* ---------------------- output of data -------------------------------- */
   fp2=fopen("c:data4.txt","w");
   fprintf(fp2,"--- input data ---\n");
   fprintf(fp2,"elas=%14.5e   pois=%14.5e\n",elas,pois);
   fprintf(fp2,"np=%5d  ne=%5d\n",np,ne);
   for(i=0;i<np;i++){
     fprintf(fp2,"%3d x[%3d]=%10.5lf y[%3d]=%10.5lf\n",i+1,i+1,x[i],i+1,y[i]);
   }
   for(i=0;i<ne;i++){
     fprintf(fp2,"%3d i=%4d j=%4d k=%4d\n",i+1,ielm[i][0],ielm[i][1],ielm[i][2]);
   }
   for(i=0;i<np;i++){
     fprintf(fp2,"node=%3d code=%3d\n",i+1,koteix[i],koteiy[i]);
   }
   for(i=0;i<np;i++){
     fprintf(fp2,"%3d loadx=%10.3e loady=%10.3e\n",i+1,loadx[i],loady[i]);
   }
   fprintf(fp2,"----- results -----\n");
   fprintf(fp2,"----- displacement of node -----\n");
   for(i=0;i<np;i++){
     fprintf(fp2,"%3d %14.5e %14.5e\n",i+1,disp[(i+1)*2-2],disp[(i+1)*2-1]);
   }
   fprintf(fp2,"----- strain of element -----\n");
   for(i=0;i<ne;i++){
     fprintf(fp2,"%3d %14.5e %14.5e %14.5e\n",i+1,eps[i][0],eps[i][1],eps[i][2]);
   }
   fprintf(fp2,"----- stress of element -----\n");
   for(i=0;i<ne;i++){
     fprintf(fp2,"%3d %14.5e %14.5e %14.5e\n",i+1,sigma[i][0],sigma[i][1],sigma[i][2]);
   }
   fclose(fp2);
   return(0);
}