/*
What's new:
C44 - [Martin Zoubek] first public release
C49 - [Martin Zoubek] * updated functions Mesh_Read, Mesh_Write and
      Mesh_FindNeighbours[1] for periodic boundary conditions
      * added missing parentheses around Mesh_EdgeOpposite
      * Mesh_Read now will set rho, v, p to 1 if no .sol given
      (for optimizing meshes by AMA)
*/

#include <stdio.h>
#include <stdlib.h>
#include "Error.h"
#include "FileIO.h"
#include "Memory.h"
#include "Mesh.h"
#include "Progress.h"
#include "Types.h"


tMesh_Mesh *Mesh_Create(void){
  tMesh_Mesh *m;
  if(!Memory_GetM(m,tMesh_Mesh,1))
    return 0;
  m->Point=0;
  m->VolumeElement=m->VolumeElementLast=0;
  m->SurfaceElement=0;
  /*  m->Red=0; */
  m->npoin=m->nelem=m->nbound=m->nred=0;
  m->NeighboursFound=FALSE;
  m->Surface=0;
  return m;
  }

void Mesh_Destroy(tMesh_Mesh *m){
  if(!m)
    return;
  while(m->VolumeElement)
    Mesh_VolumeElementDelete(m,m->VolumeElement);
  while(m->SurfaceElement)
    Mesh_SurfaceElementDelete(m,m->SurfaceElement); 
  while(m->Point)
    Mesh_PointDelete(m,m->Point);
  Memory_Free(m);
  }

tMesh_Point *Mesh_PointAdd(tMesh_Mesh *m,double x,double y){
  tMesh_Point *p;
  if(!m||!Memory_GetM(p,tMesh_Point,1))
    return 0;
  p->Point[0]=x;
  p->Point[1]=y;
  p->Ball=0;
  p->SurfaceBall=0;
  p->Prev=0;
  p->Hanging=FALSE;
  p->PeriodicCopy=0;
  p->Next=m->Point;
  if(p->Next)
    p->Next->Prev=p;
  m->Point=p;
  m->npoin++;
  return p;
  }

void Mesh_PointDelete(tMesh_Mesh *m,tMesh_Point *p){
  if(!m||!p)
    return;
  if(p->Next)
    p->Next->Prev=p->Prev;
  if(p->Prev)
    p->Prev->Next=p->Next;
  else
    m->Point=p->Next;
  Memory_Free(p);
  m->npoin--;
  }

tMesh_VolumeElement *Mesh_VolumeElementAdd(tMesh_Mesh *m,tMesh_Point *p1,
  tMesh_Point *p2,tMesh_Point *p3,
  double w11,double  w12,double w13,double  w14,double w21,double  w22,
  double w23,double  w24,double w31,double  w32,double w33,double  w34){
  tMesh_VolumeElement *ve;
  uint i,j;
  #define Swap(p1,p2){\
    tMesh_Point *_p;\
    _p=p1;\
    p1=p2;\
    p2=_p;\
    }
  if(!m||!Memory_GetM(ve,tMesh_VolumeElement,1))
    return 0;

  /*  if((pointer)p1>(pointer)p2)  //sort values that pi<pj for all i<j
    Swap(p1,p2);
  if((pointer)p1>(pointer)p3)
    Swap(p1,p3);
  if((pointer)p2>(pointer)p3)
    Swap(p2,p3); */

  ve->P[0]=p1;
  ve->P[1]=p2;
  ve->P[2]=p3;
  ve->w[0][0] = w11;
  ve->w[0][1] = w12;
  ve->w[0][2] = w13;
  ve->w[0][3] = w14;
  ve->w[1][0] = w21;
  ve->w[1][1] = w22;
  ve->w[1][2] = w23;
  ve->w[1][3] = w24;
  ve->w[2][0] = w31;
  ve->w[2][1] = w32;
  ve->w[2][2] = w33;
  ve->w[2][3] = w34;

  /*        printf("Init: %8.5g %8.5g %8.5g %8.5g %8.5g %8.5g %8.5g %8.5g \n",
	    w11,w12,w12,w14,w21,w22,w23,w24); */

  /*  printf("Nodes=%i,%i,%i",ve->P[0],ve->P[1],ve->P[2]);
  printf("Density, pressure=%0.10g,%0.10g\n",ve->rho,ve->p); */


  for(i=0;i<3;i++){  //enqueue in balls
    ve->BallPrev[i]=0;
    ve->BallNext[i]=ve->P[i]->Ball;
    if(ve->BallNext[i]){
      for(j=0;ve->BallNext[i]->P[j]!=ve->P[i];j++);
      ve->BallNext[i]->BallPrev[j]=ve;
      }
    ve->P[i]->Ball=ve;
    }
  ve->Last=FALSE;
  ve->Prev=0;
  ve->Next=m->VolumeElement;
  if(ve->Next)
    ve->Next->Prev=ve;
  else
    m->VolumeElementLast=ve;
  m->VolumeElement=ve;
  m->nelem++;
  return ve;
  #undef Swap
  }

void Mesh_VolumeElementUnqueue(tMesh_Mesh *m,tMesh_VolumeElement *ve){
  uint i,j;
  if(!m||!ve)
    return;
  for(i=0;i<4;i++){  //remove element from the ball queues
    if(ve->BallNext[i]){
      for(j=0;ve->BallNext[i]->P[j]!=ve->P[i];j++);
      ve->BallNext[i]->BallPrev[j]=ve->BallPrev[i];
      }
    if(ve->BallPrev[i]){
      for(j=0;ve->BallPrev[i]->P[j]!=ve->P[i];j++);
      ve->BallPrev[i]->BallNext[j]=ve->BallNext[i];
      }
    else
      ve->P[i]->Ball=ve->BallNext[i];
    }
  if(ve->Next)  //remove from the main queue
    ve->Next->Prev=ve->Prev;
  else
    m->VolumeElementLast=ve->Prev;
  if(ve->Prev)
    ve->Prev->Next=ve->Next;
  else
    m->VolumeElement=ve->Next;
  m->nelem--;
  }

void Mesh_VolumeElementDelete(tMesh_Mesh *m,tMesh_VolumeElement *ve){
  if(!m||!ve)
    return;
  Mesh_VolumeElementUnqueue(m,ve);
  Memory_Free(ve);
  }

void Mesh_VolumeElementLast(tMesh_Mesh *m,tMesh_VolumeElement *ve){
  if(!m||!ve)
    return;
  ve->Last=TRUE;
  if(!ve->Next)
    return;
  if(ve->Prev)
    ve->Prev->Next=ve->Next;
  else
    m->VolumeElement=ve->Next;
  ve->Next->Prev=ve->Prev;
  m->VolumeElementLast->Next=ve;
  ve->Prev=m->VolumeElementLast;
  ve->Next=0;
  m->VolumeElementLast=ve;
  }

tMesh_SurfaceElement *Mesh_SurfaceElementAdd(tMesh_Mesh *m,tMesh_Point *p1,
  tMesh_Point *p2,uint surfnr){
  tMesh_SurfaceElement *se;
  uint i,j;
  #define Swap(p1,p2){\
    tMesh_Point *_p;\
    _p=p1;\
    p1=p2;\
    p2=_p;\
    }
  if(!m||!Memory_GetM(se,tMesh_SurfaceElement,1))
    return 0;
  /*  if((pointer)p1>(pointer)p2)  //sort points
    Swap(p1,p2); */
  se->P[0]=p1;
  se->P[1]=p2;
  se->surfnr=surfnr;
  for(i=0;i<2;i++){  //enqueue in balls
    se->SurfaceBallPrev[i]=0;
    se->SurfaceBallNext[i]=se->P[i]->SurfaceBall;
    if(se->SurfaceBallNext[i]){
      for(j=0;se->SurfaceBallNext[i]->P[j]!=se->P[i];j++);
      se->SurfaceBallNext[i]->SurfaceBallPrev[j]=se;
      }
    se->P[i]->SurfaceBall=se;
    }
  se->Prev=0;
  se->Next=m->SurfaceElement;
  if(se->Next)
    se->Next->Prev=se;
  m->SurfaceElement=se;
  m->nbound++;
  return se;
  #undef Swap
  }

void Mesh_SurfaceElementDelete(tMesh_Mesh *m,tMesh_SurfaceElement *se){
  uint i,j;
  if(!m||!se)
    return;
  for(i=0;i<3;i++){  //remove element from the ball queues
    if(se->SurfaceBallNext[i]){
      for(j=0;se->SurfaceBallNext[i]->P[j]!=se->P[i];j++);
      se->SurfaceBallNext[i]->SurfaceBallPrev[j]=se->SurfaceBallPrev[i];
      }
    if(se->SurfaceBallPrev[i]){
      for(j=0;se->SurfaceBallPrev[i]->P[j]!=se->P[i];j++);
      se->SurfaceBallPrev[i]->SurfaceBallNext[j]=se->SurfaceBallNext[i];
      }
    else
      se->P[i]->SurfaceBall=se->SurfaceBallNext[i];
    }
  if(se->Next)
    se->Next->Prev=se->Prev;
  if(se->Prev)
    se->Prev->Next=se->Next;
  else
    m->SurfaceElement=se->Next;
  Memory_Free(se);
  m->nbound--;
  }

/* tMesh_Red *Mesh_RedAdd(tMesh_Mesh *m,tMesh_VolumeElement *ve){
  tMesh_Red *r;
  if(!m||!Memory_GetM(r,tMesh_Red,1))
    return 0;
  r->Element=ve;
  r->Prev=0;
  r->Next=m->Red;
  if(r->Next)
    r->Next->Prev=r;
  m->Red=r;
  m->nred++;
  return r;
  } */

/* void Mesh_RedDelete(tMesh_Mesh *m,tMesh_Red *r){
  if(!m||!r)
    return;
  if(r->Next)
    r->Next->Prev=r->Prev;
  if(r->Prev)
    r->Prev->Next=r->Next;
  else
    m->Red=r->Next;
  Memory_Free(r);
  m->nred--;
  }  */

/* read mesh from file(s), .sol file must not exist (sol is 0) */
tMesh_Mesh *Mesh_Read(FILE *grid,FILE *sol){
  tMesh_Mesh *m;
  uint npoin,nelem,nbound,nper,counter=0;
  tMesh_Point **point,*p1,*p2;
  uint i,j,k;
  float r;
  tMesh_VolumeElement *ve;
 tMesh_SurfaceElement *se;
  if(!grid||(m=Mesh_Create())==0)
    return 0;
  npoin=FileIO_ReadInteger(grid);
  if(!Memory_GetM(point,tMesh_Point *,npoin))
    goto b02;
  nelem=FileIO_ReadInteger(grid);
  nbound=FileIO_ReadInteger(grid);

  /* pomocne nacitani, nbc */
  i=FileIO_ReadInteger(grid);
  /* konec pomocneho nacitani */
  /*  nacitani periodicity */
  m->nper=0;
  for(j=0;j<2;j++){
    for(k=0;k<2;k++)
      m->rper[2*j+k]=FileIO_ReadFloat(grid);
    for(k=0;k<2;k++)
      m->iper[2*j+k]=FileIO_ReadInteger(grid);
  }
  for(j=0;j<2;j++){
    if(m->rper[2*j] != 0. || m->rper[2*j+1] != 0.){
      if(m->iper[2*j] == 0 && m->iper[2*j+1] == 0){
	printf("Mischmatch in periodicity parametres %7.5g %7.5g %d %d\n",
	       m->rper[2*j],m->rper[2*j+1], m->iper[2*j],m->iper[2*j+1]);
	exit;
      }
      m->nper++;
    }
  }
  for(j=0;j<m->nper;j++){
    printf("Periodicity params %d: (%7.5g,%7.5g) vector, %d %d (boundaries)\n",
	   m->nper,m->rper[2*j],m->rper[2*j+1], m->iper[2*j],m->iper[2*j+1]);
  }

  printf("%i points, %i volume element%s, %i surface elements\n",npoin,
    nelem,nelem>1?"s":"",nbound);
  /*  if(sol)
    m->kappa=FileIO_ReadFloat(sol);
  else 
    m->kappa=1.4; */

  Progress_Init("Reading mesh",npoin+nelem+nbound);
  for(i=0;i<npoin;i++){  //read points
    if((point[i]=Mesh_PointAdd(m,FileIO_ReadFloat(grid),
        FileIO_ReadFloat(grid)))==0)
      goto b01;
    Progress_Update(++counter);
    }

  for(i=0;i<nelem;i++){  //read elements
    int k1,k2,k3;
    k1 = FileIO_ReadInteger(grid)-1;
    k2 = FileIO_ReadInteger(grid)-1;
    k3 = FileIO_ReadInteger(grid)-1;
    //  NEW TRANSLATOR doesn't compute ????
    if(sol){
      /*
      ve=Mesh_VolumeElementAdd(m,point[FileIO_ReadInteger(grid)-1],
        point[FileIO_ReadInteger(grid)-1],point[FileIO_ReadInteger(grid)-1],
        FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),
        FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),
        FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),
        FileIO_ReadFloat(sol),FileIO_ReadFloat(sol),FileIO_ReadFloat(sol));
      */
      double r[12];
      int j;
      for(j=0; j<12; j++)
	r[j]=FileIO_ReadFloat(sol);
      ve=Mesh_VolumeElementAdd(m,point[k1],point[k2],point[k3],
                               r[0],r[1],r[2],r[3],r[4],r[5],r[6],r[7],r[8],
			       r[9],r[10],r[11]);
      
      }
    else{
    //  NEW TRANSLATOR doesn't compute ????
      /*
      ve=Mesh_VolumeElementAdd(m,point[FileIO_ReadInteger(grid)-1],
        point[FileIO_ReadInteger(grid)-1],point[FileIO_ReadInteger(grid)-1],
        1,1,1,1,1,1,1,1,1,1,1,1);
      */
      ve=Mesh_VolumeElementAdd(m,point[k1],point[k2],point[k3],
			       1,1,1,1,1,1,1,1,1,1,1,1);
	
    }
    
    if(Triangle_Deter(ve->P[0]->Point,ve->P[1]->Point,ve->P[2]->Point) <= 0){
      printf("%7.5g %7.5g %7.5g %7.5g  %7.5g %7.5g elements ZERO %7.5g  \n",
	     ve->P[0]->Point[0],ve->P[0]->Point[1],
	     ve->P[1]->Point[0],ve->P[1]->Point[1],
	     ve->P[2]->Point[0],ve->P[2]->Point[1],
	     Triangle_Deter(ve->P[0]->Point,ve->P[1]->Point,ve->P[2]->Point)); 
    } 

    if(!ve)
      goto b01;
    Progress_Update(++counter);
  }


  /*  ve=m->VolumeElement; 
  i=0;
  while(ve){
    i++;
        printf("Elementy: %i %10.5g \n", i, ve->v[0]);  

    ve=ve->Next;
    }
*/



  for(i=0;i<nbound;i++){  //read surface elements
    int k1,k2,k3;
    k1 = FileIO_ReadInteger(grid)-1;
    k2 = FileIO_ReadInteger(grid)-1;
    k3 = FileIO_ReadInteger(grid);
    //  NEW TRANSLATOR doesn't compute ????
    /*
    if(Mesh_SurfaceElementAdd(m,point[FileIO_ReadInteger(grid)-1],
           point[FileIO_ReadInteger(grid)-1],FileIO_ReadInteger(grid))==0)
    */
    if(Mesh_SurfaceElementAdd(m,point[k1],point[k2],k3)==0)
      goto b01;
    Progress_Update(++counter);
    }
  Progress_Done();

  /*
  se=m->SurfaceElement;
  while(se){
    printf("%7.5g %7.5g %7.5g %7.5g %d  %d  %d Curved \n",
	   se->P[0]->Point[0],se->P[0]->Point[1],
	   se->P[1]->Point[0],se->P[1]->Point[1],se->surfnr,se->curveelement,
	   nbound);
    
    se=se->Next;
  }
  */

  nper=FileIO_ReadInteger(grid);
  if(nper){
    printf("%i periodic point pair%s\n",nper,nper>1?"s":"");
    Progress_Init("Reading periodic point pairs",nper);
    counter=0;
    for(i=0;i<nper;i++){  //read periodic point pairs
      p1=point[FileIO_ReadInteger(grid)-1];
      p2=point[FileIO_ReadInteger(grid)-1];
      p1->PeriodicCopy=p2;
      p2->PeriodicCopy=p1;
      Progress_Update(++counter); 
      /*      printf("%7.5g %7.5g\n",p1->Point[0],p1->Point[1]);
	      printf("%7.5g %7.5g\n\n",p2->Point[0],p2->Point[1]);
      */
      }
    Progress_Done();
    }
  Memory_Free(point);
  return m;
b01:
  Memory_Free(point);
b02:
  Mesh_Destroy(m);
  return 0;
  }

  tMesh_Point **vMesh_Write_QuicksortBase;

static int Mesh_Write_QuicksortCompareFunction(const void *a,const void *b){
  Progress_Update((tMesh_Point **)a-vMesh_Write_QuicksortBase);
  return (pointer)*(tMesh_Point **)a-(pointer)*(tMesh_Point **)b;
  }

//write mesh, .grid xor .sol must not exist
void Mesh_Write(tMesh_Mesh *m,FILE *grid,FILE *sol){
  uint i,j,counter=0,p1,p2,p3,nper=0;
  tMesh_Point *p,**point;
  tMesh_VolumeElement *ve;
  tMesh_SurfaceElement *se;
  #define FindPoint(p,i){\
    uint _a,_b;\
    _a=0;\
    _b=m->npoin-1;\
    while(point[i=(_a+_b)>>1]!=(p))\
      if((pointer)point[i]<(pointer)(p))\
        _a=i+1;\
      else\
        _b=i-1;\
    }
  if(!m||(!grid&&!sol))
    return;
  if(!grid){  //warning: the mesh must not be changed!!!
       /*    fprintf(sol,"%i\n%0.10g\n",m->nelem,m->kappa); */
    Progress_Init("Writing mesh",m->nelem);
    ve=m->VolumeElement;
    while(ve->Next)
      ve=ve->Next;
    while(ve){
      for(j=0;j<3;j++)
        fprintf(sol,"%0.16g %0.16g %0.16g  %0.16g\n",
          ve->w[j][0],ve->w[j][1],ve->w[j][2],ve->w[j][3]);

      /*      fprintf(sol,"%0.10g %0.10g %0.10g  %0.10g\n",ve->rho,ve->v[0],
        ve->v[1],ve->p); */ 
      Progress_Update(++counter);
      ve=ve->Prev;
      }
    Progress_Done();
    return;
    }
  if(!Memory_GetM(point,tMesh_Point *,m->npoin))  //sort points
    return;
  p=m->Point;
  i=0;
  while(p){
    point[i++]=p;
    if(p->PeriodicCopy)
      nper++;
    p=p->Next;
    }
  nper>>=1;
  Progress_Init("Sorting points",m->npoin);
  vMesh_Write_QuicksortBase=point;
  qsort((void *)point,m->npoin,sizeof(tMesh_Point *),
    Mesh_Write_QuicksortCompareFunction);
  Progress_Done();
  Progress_Init("Writing mesh",m->npoin+m->nelem+m->nbound+nper);
  fprintf(grid,"%i\n%i\n%i\n",m->npoin,m->nelem,m->nbound);
  /*  if(sol)
    fprintf(sol,"%i\n%0.10g\n",m->nelem,m->kappa); */
  for(i=0;i<m->npoin;i++){
    fprintf(grid,"%0.10g %0.10g %0.10g\n",point[i]->Point[0],
      point[i]->Point[1],point[i]->Point[2]);
    Progress_Update(++counter);
    }
  ve=m->VolumeElement;
  while(ve){
    FindPoint(ve->P[0],p1);
    FindPoint(ve->P[1],p2);
    FindPoint(ve->P[2],p3);
    fprintf(grid,"%i %i %i %i\n",p1+1,p2+1,p3+1);
    if(sol)
      for(j=0;j<3;j++)
        fprintf(sol,"%0.10g %0.10g %0.10g  %0.10g\n",
          ve->w[j][0],ve->w[j][1],ve->w[j][2],ve->w[j][3]);

    /* fprintf(sol,"%0.10g %0.10g %0.10g  %0.10g\n",ve->rho,ve->v[0]*ve->rho,
        ve->v[1]*ve->rho, ve->p/(m->kappa-1)+
      0.5*ve->rho*double2_norm2(ve->v) ); */

    /*      fprintf(sol,"%0.10g %0.10g %0.10g %0.10g\n",ve->rho,ve->v[0],
        ve->v[1],ve->p);*/
    Progress_Update(++counter);
    ve=ve->Next;
    }
  se=m->SurfaceElement;
  while(se){
    FindPoint(se->P[0],p1);
    FindPoint(se->P[1],p2);
    fprintf(grid,"%i %i %i\n",p1+1,p2+1,se->surfnr);
    Progress_Update(++counter);
    se=se->Next;
    }
  fprintf(grid,"%i\n",nper);
  for(i=0;i<m->npoin;i++)
    if(point[i]->PeriodicCopy&&
      (pointer)(point[i])<(pointer)(point[i]->PeriodicCopy)){
      FindPoint(point[i]->PeriodicCopy,p1);
      fprintf(grid,"%i %i\n",i+1,p1+1);
      Progress_Update(++counter);
      }
  Progress_Done();
  Memory_Free(point);
  }

void Mesh_FindNeighbours1(tMesh_VolumeElement *ve){
  uint i,j;
  tMesh_VolumeElement *ve2;
  tMesh_SurfaceElement *se;
  static uint perm[3][2]={{0,1},{1,2},{2,0}};
  for(i=0;i<3;i++){
    ve->Neighbour[i]=0;
    ve2=ve->P[perm[i][0]]->Ball;
    while(ve2){
      if(ve!=ve2&&Mesh_ContainsEdge(ve2,ve->P[perm[i][1]],
        ve->P[perm[i][2]])){
        ve->Neighbour[i]=ve2;
        break;
        }
      for(j=0;ve2->P[j]!=ve->P[perm[i][0]];j++);
      ve2=ve2->BallNext[j];
      }
    if(!ve->Neighbour[i]&&ve->P[perm[i][0]]->PeriodicCopy&&
      ve->P[perm[i][1]]->PeriodicCopy&&ve->P[perm[i][2]]->PeriodicCopy){
      ve2=ve->P[perm[i][0]]->PeriodicCopy->Ball;
      while(ve2){
        if(Mesh_ContainsEdge(ve2,ve->P[perm[i][1]]->PeriodicCopy,
          ve->P[perm[i][2]]->PeriodicCopy)){
          ve->Neighbour[i]=ve2;
          break;
          }
        for(j=0;ve2->P[j]!=ve->P[perm[i][0]]->PeriodicCopy;j++);
        ve2=ve2->BallNext[j];
        }
      }
    ve->Surface[i]=0;
    se=ve->P[perm[i][0]]->SurfaceBall;
    while(se){
      if(Mesh_SurfaceContainsEdge(se,ve->P[perm[i][1]],ve->P[perm[i][2]])){
        ve->Surface[i]=se;
        break;
        }
      for(j=0;se->P[j]!=ve->P[perm[i][0]];j++);
      se=se->SurfaceBallNext[j];
      }
    }
  }

void Mesh_FindNeighbours(tMesh_Mesh *m){
  tMesh_VolumeElement *ve;
  uint counter=0;
  if(!m)
    return;
  if(m->NeighboursFound)
    return;
  Progress_Init("Looking up neighbours",m->nelem);
  ve=m->VolumeElement;
  while(ve){
    ve->Parent=0;
    Mesh_FindNeighbours1(ve);
    Progress_Update(++counter);
    ve=ve->Next;
    }
  Progress_Done();
  m->NeighboursFound=TRUE;
  }

