#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <string.h>
#include <math.h>


#define MAXMAT 1000
#define MAXCND 1000
   
char projname[80];
int i, ielem, inod, icnd;
double *x, *y;
int *N, *imat;
int  nodc[MAXCND];
double rho[MAXMAT], wval[MAXCND];
int Nelem, Nnod, Nmat, Ncnd;
double x_CG, y_CG;


void input(void);
void calculate(void);
void output(void);

void main (int argc, char *argv[])       
 {      

   strcpy (projname, argv[1]);
   input();
   calculate();
   output();
 }

void input ()
 {
   char filedat[80], fileerr[80], sau1[80], sau2[80];
   FILE *fp, *ferr;
   int aux,j, error=0;

   void jumpline (FILE*);  
   

   strcpy(filedat, projname);
   strcat(filedat,".dat");
   fp = fopen(filedat, "r");
   


   for (i=0;  i<6; i++)
     jumpline (fp);
   fscanf(fp, "%d %d", &Nelem, &Nnod);
 
    
   x=(double *) malloc((Nnod+1)*sizeof(double)); if (x==NULL) {error=1;}
   y=(double *) malloc((Nnod+1)*sizeof(double)); if (y==NULL) {error=1;}

   N=   (int *) malloc((Nelem+1)*3*sizeof(int)); if (N==NULL) {error=1;}
   imat=(int *) malloc((Nelem+1)*sizeof(int));   if (N==NULL) {error=1;}
   
   if (error)
   {
      strcpy(fileerr, projname);
      strcat(fileerr,".err");
      ferr = fopen(fileerr, "w");
      fprintf(ferr, "\n \n \n *****  ERROR: Not enough memory. ***** \n");               
      fprintf(ferr, "(Try to calculate with less elements)\n");               

      exit(1);  
    }

   for (i=0;  i<6; i++)
     jumpline (fp);

   /* lectura coordenadas nodos */
   for (inod=1; inod<=Nnod; inod++)       
       fscanf (fp, "%d %lf %lf", &aux, &x[inod], &y[inod]);  



   for (i=0;  i<6; i++)
     jumpline (fp);
  
/* lectura conectividades elementos */

   for (ielem=1; ielem<=Nelem; ielem++){
     fscanf (fp, "%d", &aux);
     for(j=0;j<3;j++) fscanf (fp, "%d", &N[(ielem-1)*3+j]);
     fscanf (fp, "%d", &imat[ielem]);
     if (imat[ielem]==0)    
         {
          strcpy(fileerr, projname);
          strcat(fileerr,".err");
          ferr = fopen(fileerr, "w");
          fprintf(ferr, "\n \n \n **ERROR: Elements with no material!!** \n"); 
          exit(1);  
         }
   }


   for (i=0;  i<5; i++)
     jumpline (fp);
   fscanf(fp, "%s %s  %d",sau1, sau2, &Nmat );

   for (i=0;  i<3; i++)
     jumpline (fp);

   /* lectura densidad materiales */
   for (i=1; i<=Nmat; i++)
       fscanf (fp, "%d %lf", &aux, &rho[i]);
       
    
   for (i=0;  i<4; i++)
     jumpline (fp);

   /* lectura condiciones puntuales de peso */
   fscanf(fp, "%d", &Ncnd);  

   for (i=0;  i<6; i++)
     jumpline (fp);

   for (icnd=1; icnd<=Ncnd; icnd++)
     { 
       fscanf (fp, "%d %lf", &nodc[icnd], &wval[icnd]);
       jumpline (fp);
       
       
     }

   fclose (fp);
 }

void calculate ()
 {
   double v;
   int n1, n2, n3;
   int mat;
   double x_CGi, y_CGi;
   double x_num=0, y_num=0,  den=0;


   
   /* para cada elemento del mallado... */
   for (ielem=1; ielem<=Nelem; ielem++)
     {
       n1= N[0+(ielem-1)*3];
       n2= N[1+(ielem-1)*3];
       n3= N[2+(ielem-1)*3];

       /* calculo el volumen (caso 2D coincide con el area) */
       
      v=fabs(x[n1]*y[n2]+x[n2]*y[n3]+x[n3]*y[n1]-x[n1]*y[n3]-x[n2]*y[n1]-x[n3]*y[n2])/2;

   

       /* calculo el centro geometrico del elemento
	  (que coincide con el centro de gravedad del elemento siendo
	  la densidad uniforme en su interior) */
       x_CGi= (x[n1]+x[n2]+x[n3])/3;
       y_CGi= (y[n1]+y[n2]+y[n3])/3;

    /* actualizo los sumatorios */
       mat= imat[ielem];
       x_num+= rho[mat]*v*x_CGi;
       y_num+= rho[mat]*v*y_CGi;
 
       den+= rho[mat]*v;
     }

   /* influencia pesos (condiciones) puntuales */
   for (icnd=1; icnd<=Ncnd; icnd++)
     {
       inod= nodc[icnd];
       x_num+= wval[icnd]*x[inod];
       y_num+= wval[icnd]*y[inod];
       
       den+= wval[icnd]; 
     }  

  
   
   x_CG= (x_num/den);
   y_CG= (y_num/den);
 
      
 }

void output()
 {
   char filedat[80];
   FILE *fp, *fptest;
   double raiz;
   double pot;

   strcpy(filedat, projname);
   strcat(filedat,".boh");
   fptest = fopen(filedat, "w");
   fprintf (fptest, "FILE: %s\n", projname);
   fprintf (fptest, "CMAS2D\n2D routine to calculate the mass center of an heterogeneus object.\nJanuary 2000\tRienzi Gianfranco\t\t\t\tCIMNE\n");
   fprintf (fptest, "\n\n\t====> mass center: %lf %lf \n", x_CG, y_CG);

    /* escritura en .flavia.dat */
   strcpy(filedat, projname);
   strcat(filedat,".flavia.dat");
   fp = fopen(filedat, "w");
   fprintf (fp, "FILE: %s\n", projname);
   fprintf (fp, "CMAS\n");
   fprintf (fp, "Programa para el calculo del centro de masas\n");
   fprintf (fp, "\t\t\t\t\tCIMNE\n");
   fprintf (fp, "2000 G. Rienzi\n");
   fprintf (fp, "\n");
  
   fprintf (fp, " %6d %6d      3\n", Nelem, Nnod);   

   fprintf (fp, "\t\tNODES\n");   
   for (inod=1; inod<=Nnod; inod++)       
     fprintf (fp, " %6d %14.6e %14.6e\n", inod, x[inod], y[inod]);  
   fprintf (fp, "\t\tCONECTIVITIES\n");   
   
   
    for (ielem=1; ielem<=Nelem; ielem++){
     fprintf (fp, " %6d %6d %6d %6d %6d\n", ielem, N[0+(ielem-1)*3], N[1+(ielem-1)*3], N[2+(ielem-1)*3], imat[ielem]);  
   }
   
   
   fclose (fp);
  
   /* escritura en .flavia.res */
   strcpy(filedat, projname);
   strcat(filedat,".flavia.res");
   fp = fopen(filedat, "w");
   fprintf (fp, "MC-DISTANCE    2    1    1    1    0\n");
   for (inod=1; inod<=Nnod; inod++)
     {
       /* calculo la distancia de cada nodo dal centro de masa ("potencial") */
       raiz= (x_CG-x[inod])*(x_CG-x[inod]) + (y_CG-y[inod])*(y_CG-y[inod]);
       pot = sqrt (raiz);
       fprintf (fp, " %6d %14.6lf\n", inod, pot);
     }
   fclose (fp);   
   fclose (fptest);
   free(x);
   free(y);
   free(N);
   free(imat);

 }

/* funcion de salto de linea */
void jumpline (FILE* filep)
 {
   char buffer[80];
   
   fgets(buffer, 80, filep);
 }

