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

typedef double Daten_Typ;
typedef unsigned int Index_Typ;
typedef unsigned int Zaehler_Typ;
typedef Zaehler_Typ Anz_Typ;

#define BSize 1024 
char NOut[BSize],NIn[BSize];
Daten_Typ Global_Time=0.0;

#ifdef PART3D
#define DIM 3 
#define NpC 8
#define CpM 8  
const char CellStr[]="hex 8"; 
#else
#define DIM 2 
#define NpC 4 
#define CpM 4 
const char CellStr[]="quad 4";
#endif

#ifdef BOUSS
#define Num_Field 3
const char * Field_Names[Num_Field]=
{
  "pressure 1",
  "streamfunction 1",
  "temperature 1"
};

const char * Tracer_Names[Num_Field]=
{
  "Druck",
  "Stromfkt",
  "Temp"
};
#else
#define Num_Field 2 
const char *Field_Names[Num_Field]=
{
  "pressure 1",
  "streamfunction 1"
};

const char *Tracer_Names[Num_Field]=
{
  "Druck",
  "Stromfkt"
};
#endif

typedef struct
{
  Anz_Typ Anz_Index;
  Index_Typ *I;
} Index_List;

typedef struct
{
  Anz_Typ Anz_Nodes;
  Daten_Typ *K[DIM];
} Nodes_Typ;

/* Moegliche Zustaende fuer Reg_Flag[i] */
#define REG_UNKNOWN 0 
#define REG_REGULAR 1 
#define REG_IRREGULAR 2 

typedef struct
{
  Anz_Typ Anz_Cells;
  Index_Typ *Vert;
  unsigned char *Reg_Flag;
  /* Formfaktoren fuer verbesserte Interpolation */
  /* Daten_Typ *Form; */
} Cells_Typ;

typedef struct 
{
  Nodes_Typ Nodes;
  Cells_Typ Cells;	
} Gitter_Typ;

typedef struct
{
  Gitter_Typ Gitter;
  Daten_Typ *V[DIM];
  Daten_Typ *Field[Num_Field];
} Input_Typ;

typedef struct
{
  Daten_Typ L[NpC];
} B_Coord_Typ;

typedef struct
{
  Anz_Typ Anz_Tracer;
  Daten_Typ *K[DIM];
  Daten_Typ *Farbe,*Expire;
  Index_Typ *InCell;
  B_Coord_Typ *B_Coords;
} Tracer_Typ;

typedef struct
{
  Gitter_Typ Gitter;
  Tracer_Typ Tracers;
  Daten_Typ *Field[Num_Field];
} Output_Typ;

typedef struct
{
  Daten_Typ P[2][DIM],Farbe,Haltbarkeit;
  Anz_Typ N[DIM];
  Anz_Typ Modul,Offset,Start,Stop;
} Block_Typ;

Output_Typ O;
Input_Typ I;
/* Die Einzelnen Ebenen des Mehrgitters */
Anz_Typ Anz_Level=0;
Gitter_Typ *MG;
/* Liste der "entarteten" Zellen */
Index_List *CL;

const Daten_Typ EPS=0.5E-7;/* 0.2E-4; */
const Daten_Typ Small_EPS=0.5E-7;
 
void Lese_Daten(FILE *f,Anz_Typ Anz,Daten_Typ *Daten)
{
  Zaehler_Typ i;
  if(Daten==NULL)
    {
      for(i=0;i<Anz;i++) fscanf(f,"%*f");
    }
  else
    {
      for(i=0;i<Anz;i++,Daten++) fscanf(f,"%lf",Daten);
    }
}
void Schreibe_Daten(FILE *f,Anz_Typ Anz,Daten_Typ *Daten)
{
  Zaehler_Typ i;
  if(Daten==NULL)
    {
      for(i=0;i<Anz;i++) fprintf(f,"%.8f\n",0.0);
    }
  else
    {
      for(i=0;i<Anz;i++,Daten++) fprintf(f,"%.8f\n",*Daten);
    }
}
void Lese_Index(FILE *f,Anz_Typ Anz,Index_Typ *Index)
{
  Zaehler_Typ i;
  if(Index==NULL)
    {
      for(i=0;i<Anz;i++) fscanf(f,"%*u");
    }
  else
    {
      for(i=0;i<Anz;i++,Index++) fscanf(f,"%u",Index);
    }
}
void Schreibe_Index(FILE *f,Anz_Typ Anz,Index_Typ *Index)
{
  Zaehler_Typ i;
  if(Index==NULL)
    {
      for(i=0;i<Anz;i++) fprintf(f,"%u\n",0);
    }
  else
    {
      for(i=0;i<Anz;i++,Index++) fprintf(f,"%u\n",*Index);
    }
}
void Lese_Anz(FILE *f,Anz_Typ *Anz)
{
  if(Anz==NULL) fscanf(f,"%*u");
  else fscanf(f,"%u",Anz);
}
void Schreibe_Anz(FILE *f,Anz_Typ Anz)
{
  fprintf(f,"%u\n",Anz);
}

void Scan_Name(FILE *f,const char *S)
{
  fscanf(f,"\n");
  fscanf(f,S);
}
 
void Schreibe_Gitter(FILE *f,Gitter_Typ *Gitter)
{
  Zaehler_Typ i;
  Anz_Typ A;
  fprintf(f,"gmvinput ascii\nnodes ");
  A=Gitter->Nodes.Anz_Nodes;
  Schreibe_Anz(f,A);
  for(i=0;i<3;i++)
    {
      if(i<DIM)Schreibe_Daten(f,A,Gitter->Nodes.K[i]+1);
      else Schreibe_Daten(f,A,NULL);
    }
  fprintf(f,"cells ");
  A=Gitter->Cells.Anz_Cells;
  Schreibe_Anz(f,A);
  for(i=1;i<=A;i++)
  { 
    fprintf(f,"%s\n",CellStr);
    Schreibe_Index(f,NpC,Gitter->Cells.Vert+i*NpC);
  }
}

void Lese_Gitter(FILE *f,Gitter_Typ *Gitter)
{
   Anz_Typ A;
   Zaehler_Typ i;
   fscanf(f,"gmvinput ascii");
   Scan_Name(f,"nodes");
   Lese_Anz(f,&A);
   if(Gitter!=NULL)
   {
     printf("Ich alloziere nun Speicher für %u Knoten...\n",A);
     Gitter->Nodes.Anz_Nodes=A;
     for(i=0;i<DIM;i++) 
       Gitter->Nodes.K[i]=
	 (Daten_Typ *)calloc(A+1,sizeof(Daten_Typ));
     printf("Ich lese nun %u Knoten ein...\n",A);
     for(i=0;i<3;i++)
       {
	 if(i<DIM) Lese_Daten(f,A,Gitter->Nodes.K[i]+1);
	 else Lese_Daten(f,A,NULL);
       }
   }
   else
   {
     printf("Ich überlese nun %u Knoten...\n",A);
     Lese_Daten(f,3*A,NULL);
   }
   Scan_Name(f,"cells");
   Lese_Anz(f,&A);
   if(Gitter!=NULL)
   {
     printf("Ich alloziere nun Speicher für %u Zellen...\n",A);
     Gitter->Cells.Anz_Cells=A;
     Gitter->Cells.Vert=(Index_Typ *)calloc((A+1)*NpC,sizeof(Index_Typ));
     Gitter->Cells.Reg_Flag=(unsigned char *)calloc(A+1,sizeof(unsigned char));
     printf("Ich lese nun %u Zellen ein...\n",A);
     for(i=1;i<=A;i++)
     {
        Scan_Name(f,CellStr);
        Lese_Index(f,NpC,Gitter->Cells.Vert+i*NpC);
     }
   }
   else
   {
     printf("Ich überlese nun %u Zellen...\n",A);
     for(i=1;i<=A;i++)
     {
       Scan_Name(f,CellStr);
       Lese_Index(f,NpC,NULL);
     }
   }
}


void Update_Input(const char *Name,Input_Typ *Input)
{
  FILE *f;
  Anz_Typ A,i;
  f=fopen(Name,"r");
  Lese_Gitter(f,NULL);
  A=Input->Gitter.Nodes.Anz_Nodes;
  printf("Ich lese nun die Datensätze ein...\n");
  Scan_Name(f,"velocity 1");
  for(i=0;i<3;i++)
    {
      if(i<DIM) Lese_Daten(f,A,Input->V[i]+1);
      else Lese_Daten(f,A,NULL);
    }
  printf("Geschwindigkeitsfeld gelesen!\n");
  Scan_Name(f,"variable");
  for(i=0;i<Num_Field;i++)
    {
      Scan_Name(f,Field_Names[i]);
      Lese_Daten(f,A,Input->Field[i]+1);
      printf("Zusatzdatenfeld %d gelesen!\n",i+1);
    }
  fclose(f);
}

void Startup_Input(const char *Name,Input_Typ *Input)
{
  Anz_Typ A,i;
  Input->Gitter=MG[Anz_Level-1];
  A=Input->Gitter.Nodes.Anz_Nodes;
  printf("Ich alloziere Speicher für diverse Datenfelder\n");
  for(i=0;i<DIM;i++) 
    Input->V[i]=(Daten_Typ *)calloc(A+1,sizeof(Daten_Typ));
  for(i=0;i<Num_Field;i++) 
    Input->Field[i]=(Daten_Typ *)calloc(A+1,sizeof(Daten_Typ));
  Update_Input(Name,Input);
}
 
void Cleanup_Nodes(Nodes_Typ *N)
{
  Zaehler_Typ i;
  N->Anz_Nodes=0;
  for(i=0;i<DIM;i++) 
    if(N->K[i]!=NULL)
    {
      free(N->K[i]);
      N->K[i]=NULL;
    }
}

void Cleanup_Cells(Cells_Typ *C)
{
  C->Anz_Cells=0;
  if(C->Vert!=NULL)
  {
    free(C->Vert);
    C->Vert=NULL;
  }
  if(C->Reg_Flag!=NULL)
  {
    free(C->Reg_Flag);
    C->Reg_Flag=NULL;
  }
}

void Cleanup_Gitter(Gitter_Typ *G)
{
  Cleanup_Nodes(&(G->Nodes));
  Cleanup_Cells(&(G->Cells));
}

void Cleanup_Input(Input_Typ *Input)
{
  Anz_Typ i;
   if(Input!=NULL)
   {
     printf("Ich gebe nun allozierten Speicher wieder frei...\n");
     for(i=0;i<DIM;i++)
       free(Input->V[i]);
     for(i=0;i<Num_Field;i++)
       free(Input->Field[i]);
     /* Gitter wird schon bei Cleanup_Suchen freigegeben! */
   }
}

Daten_Typ My_Abs(Daten_Typ Val)
{
   return (Val<0.0 ? -Val:Val);
}

#ifdef PART3D
int Test_Ebene(Gitter_Typ *Gitter,Index_Typ P1,Index_Typ P2,Index_Typ P3,
               Daten_Typ K[3])
{
   Daten_Typ Vect1[3],Vect2[3],Normal[3],Erg;
   Zaehler_Typ i;
   for(i=0;i<DIM;i++)
   {
     Vect1[i]=Gitter->Nodes.K[i][P2]-Gitter->Nodes.K[i][P1];
     Vect2[i]=Gitter->Nodes.K[i][P3]-Gitter->Nodes.K[i][P1];
   }
   for(Erg=0.0,i=0;i<DIM;i++) Erg+=Vect1[i]*Vect1[i];
   Erg=1.0/sqrt(Erg);
   for(i=0;i<DIM;i++) Vect1[i]*=Erg;
   for(Erg=0.0,i=0;i<DIM;i++) Erg+=Vect2[i]*Vect2[i];
   Erg=1.0/sqrt(Erg);
   for(i=0;i<DIM;i++) Vect2[i]*=Erg;
   /* Berechne den auf der Flaeche stehenden Normalenvektor */
   Normal[0]=Vect1[1]*Vect2[2]-Vect1[2]*Vect2[1];
   Normal[1]=-(Vect1[0]*Vect2[2]-Vect1[2]*Vect2[0]);
   Normal[2]=Vect1[0]*Vect2[1]-Vect1[1]*Vect2[0];
   for(i=0,Erg=0.0;i<3;i++) Erg+=Normal[i]*(K[i]-Gitter->Nodes.K[i][P1]);
   return (Erg>=0.0 ?1:0);
}

int Test_Zelle(Index_Typ C,Gitter_Typ *Gitter,Daten_Typ K[DIM])
{
  const Index_Typ T[6][3]=
  {
    {0,4,1},
    {1,5,2},
    {2,6,3},
    {3,7,0},
    {0,1,3},
    {4,7,5}
  };
  Zaehler_Typ i;
  Index_Typ *D;
  D=Gitter->Cells.Vert+NpC*C;
  for(i=0;i<6;i++) 
  {
    if(!Test_Ebene(Gitter,*(D+T[i][0]),*(D+T[i][1]),*(D+T[i][2]),K)) 
      return 0;
  }
  return 1;
}

unsigned char Is_Regular(Index_Typ C,Gitter_Typ *G)
{
  const Index_Typ T[8][3]=
  {
    {1,3,4},
    {0,2,5},
    {1,6,3},
    {0,2,7},
    {0,5,7},
    {1,6,4},
    {2,5,7},
    {3,6,4},
  };
  Daten_Typ Erg;
  Zaehler_Typ i,j,k;
  if(G->Cells.Reg_Flag[C]!=REG_UNKNOWN) return G->Cells.Reg_Flag[C];
  for(i=0;i<8;i++)
  {
    for(k=0;k<3;k++)
    {
      for(j=0,Erg=0.0;j<3;j++)
      {
        
        Erg+=(G->Nodes.K[j][G->Cells.Vert[NpC*C+T[i][k]]]
              -G->Nodes.K[j][G->Cells.Vert[NpC*C+i]])
             *(G->Nodes.K[j][G->Cells.Vert[NpC*C+T[i][(k+1)%3]]]
             -G->Nodes.K[j][G->Cells.Vert[NpC*C]+i]);
      }
      if(My_Abs(Erg)>=EPS)
      {
        G->Cells.Reg_Flag[C]=REG_IRREGULAR;
        return REG_IRREGULAR;
      }
    }
  }
  G->Cells.Reg_Flag[C]=REG_REGULAR;
  return REG_REGULAR;
}
#else
int Test_Linie(Gitter_Typ *Gitter,Index_Typ P1,Index_Typ P2,Daten_Typ K[2])
{
  Daten_Typ Normal[DIM],Vect[DIM],Erg;
  Zaehler_Typ i;
  for(i=0;i<DIM;i++)
    Vect[i]=K[i]-Gitter->Nodes.K[i][P1];
  Normal[0]=-(Gitter->Nodes.K[1][P2]-Gitter->Nodes.K[1][P1]);
  Normal[1]=Gitter->Nodes.K[0][P2]-Gitter->Nodes.K[0][P1];
  
  for(Erg=0.0,i=0;i<DIM;i++) Erg+=Normal[i]*Normal[i];
  Erg=1.0/sqrt(Erg);
  for(i=0;i<DIM;i++) Normal[i]*=Erg;
  
  for(Erg=0.0,i=0;i<2;i++) Erg+=Vect[i]*Normal[i];
  /**/
  if(Erg>=-EPS) return 1;
  return 0;
}

int Test_Zelle(Index_Typ C,Gitter_Typ *Gitter,Daten_Typ K[DIM])
{
  Zaehler_Typ i;
  Index_Typ *D;
  D=Gitter->Cells.Vert+NpC*C;
  for(i=0;i<NpC;i++) 
  {
    if(!Test_Linie(Gitter,*(D+i),*(D+(i+1)%NpC),K)) return 0;
  }
  return 1;
}

unsigned char Is_Regular(Index_Typ C,Gitter_Typ *G)
{
  Daten_Typ V[NpC][2],Erg;
  Zaehler_Typ i,j;
  if(G->Cells.Reg_Flag[C]!=REG_UNKNOWN) return G->Cells.Reg_Flag[C];
  for(i=0;i<NpC;i++)
  {
    for(j=0;j<DIM;j++)
      {
	V[i][j]=G->Nodes.K[j][G->Cells.Vert[NpC*C+((i+1)%NpC)]]-
	  G->Nodes.K[j][G->Cells.Vert[NpC*C+i]];
      }
  }
  for(i=0;i<NpC;i++)
  {
    Erg=V[i][0]*V[(i+1)%NpC][0]+V[i][1]*V[(i+1)%NpC][1];
    if(My_Abs(Erg)>=EPS)
    {
      G->Cells.Reg_Flag[C]=REG_IRREGULAR;
      return REG_IRREGULAR;
    }
  }
  G->Cells.Reg_Flag[C]=REG_REGULAR;
  return REG_REGULAR;
}
#endif

#ifdef MOVBC
#include "movbc.inc"
#endif

int Datei_Existiert(const char *Name)
{
  FILE *f;
  f=fopen(Name,"r");
  if(f==NULL) return 0;
  fclose(f);
  return 1;
}

void File_Error(const char *Name)
{
  fprintf(stderr,"Datei %s nicht gefunden!\n",Name);
  exit(EXIT_FAILURE);
}

/* Suche Zelle auf Mehrgitterweise */
Index_Typ Suche_Zelle(Daten_Typ K[DIM],Index_Typ Old_Cell)
{
  Index_Typ i,Akt_Level,Akt_Cell,n;
  
#ifdef MOVBC
  /* Teste zuerst, ob das Partikel im bewegten Rand ist */
  if(Test_MovBC(K)) return 0;
#endif
  Akt_Level=Anz_Level;
  Akt_Cell=Old_Cell;
  if(Akt_Cell)
  {
    /* Vergröbere zuerst das Gitter nach bedarf max. bis Level 1*/
    while((Akt_Level>1)&&!Test_Zelle(Akt_Cell,MG+Akt_Level-1,K))
    {
      Akt_Level--;
      if(Akt_Cell>MG[Akt_Level-1].Cells.Anz_Cells)
        Akt_Cell=(Akt_Cell-MG[Akt_Level-1].Cells.Anz_Cells-1)/(CpM-1)+1;
    }
    /* Vollstaendige Suche, falls erfolglos */
    if(!Test_Zelle(Akt_Cell,MG+Akt_Level-1,K)) Akt_Cell=0;
  }
  else Akt_Level=1;
  /* Mache solange vollständige Suchen in den entarteten Zellen, 
     bis ein Einstiegsmakro gefunden wurde oder feinster Level erreicht ist.
     Worst Case Fall tritt ein, wenn ein Tracer sich aus dem Gebiet
     herausbewegt!
     */
 loop1:
  while((Akt_Level<=Anz_Level)&&(!Akt_Cell))
  {
    for(i=0;i<CL[Akt_Level-1].Anz_Index;i++)
    {
      if(Test_Zelle(CL[Akt_Level-1].I[i],MG+Akt_Level-1,K))
      {
        Akt_Cell=CL[Akt_Level-1].I[i];
        break;
      }
    }
    if(!Akt_Cell) Akt_Level++;
  }
  /* Tracer befindet sich in keiner Zelle */
  if(Akt_Level>Anz_Level) return 0;
  /* Verfeinere das Gitter Schrittweise bis zum feinsten Level */
  while(Akt_Level<Anz_Level)
  {
    n=0;
    Akt_Level++;
    if(Test_Zelle(Akt_Cell,MG+Akt_Level-1,K)) n=1;
    else
      {
	for(i=1;i<CpM;i++)
	  {
	    if(Test_Zelle((Akt_Cell-1)*(CpM-1)+MG[Akt_Level-2].Cells.Anz_Cells+i,
			  MG+Akt_Level-1,K))
	      {
		Akt_Cell=(Akt_Cell-1)*(CpM-1)+MG[Akt_Level-2].Cells.Anz_Cells+i;
		n=1;
		break;
	      }
	  }
      }
    /* Wenn es nicht in einer Zelle der aktuellen Hirarchie ist, dann kann
       immer noch in einer entarteten Zelle des naechstfeineren Gitters sein.
    */
    if(!n)
      {
	Akt_Cell=0;
	Akt_Level++;
	goto loop1;
      }
  }
  if(Akt_Level!=Anz_Level) fprintf(stderr,"Schwerer Fehler beim Suchen!\n");
  return Akt_Cell;
}

void Startup_Suchen(void)
{
  Zaehler_Typ i,j,k,l;
  Index_Typ OldCell;
  Anz_Typ *A;
  Anz_Typ n=0;
  Daten_Typ K[DIM];
  FILE *f;
  Index_List L;
  L.I=NULL;L.Anz_Index=0;
  /* Zuerst Bestimme ich die Anzahl der Verfeinerungstufen des MG */
  i=0;
  do {sprintf(NIn,"#gmv/tria%d.gmv",++i);} while(Datei_Existiert(NIn));
  Anz_Level=--i;
  MG=(Gitter_Typ *)calloc(Anz_Level,sizeof(Gitter_Typ));
  A=(Anz_Typ *)calloc(Anz_Level,sizeof(Anz_Typ));
  for(i=1;i<=Anz_Level;i++)
  {
    sprintf(NIn,"#gmv/tria%d.gmv",i);
    f=fopen(NIn,"r");
    Lese_Gitter(f,MG+i-1);
    fclose(f);
    A[i-1]=MG[i-1].Nodes.Anz_Nodes;
    if(i<Anz_Level) Cleanup_Nodes(&(MG[i-1].Nodes));
  }
  for(i=0;i<Anz_Level-1;i++) 
    {
      MG[i].Nodes=MG[Anz_Level-1].Nodes;
      MG[i].Nodes.Anz_Nodes=A[i];
    }
  free(A);
  /* Debug */
  printf("Bestimme entartete Zellen...\n");
  /**/
  /* Bestimme entartete Zellen,d.h. alle Zellen, die nicht vollständig
     in einem Makro der naechstgröberen Verfeinerungstufe liegen */
  /* Debug */
  printf("Level 1\n");
  /**/
  CL=(Index_List *)calloc(Anz_Level,sizeof(Index_List));
  /* Im gröbsten Level sind alle Zellen entartet */
  CL[0].I=(Index_Typ *)calloc(MG[0].Cells.Anz_Cells,sizeof(Index_Typ));
  CL[0].Anz_Index=MG[0].Cells.Anz_Cells;
  for(i=0;i<CL[0].Anz_Index;i++) CL[0].I[i]=i+1;
  /* Bearbeite Level 2..Anz_Level */
  for(i=2;i<=Anz_Level;i++)
  {
    /* Debug */
    printf("Level %d\n",i);
    /**/
    /* Lösche Knotenliste */
    L.Anz_Index=0;
    L.I=(Index_Typ *)realloc(L.I,0);
    /* Teste alle neuen Knoten im aktuellen Level daraufhin, ob sie nicht
       innerhalb des nächstgröberen Gitters liegen mittels
       vollständiger Suche.
    */
    OldCell=0;
    for(j=MG[i-2].Nodes.Anz_Nodes+1;j<=MG[i-1].Nodes.Anz_Nodes;j++)
    {
      for(k=0;k<DIM;k++) K[k]=MG[i-1].Nodes.K[k][j];
      /*
      n=0;
      for(k=1;k<=MG[i-2].Cells.Anz_Cells;k++)
        if(Test_Zelle(k,MG+i-2,K))
        {
          n=1;
          break;
        }
	*/
      /* Experimentell */ 
      l=Anz_Level;
      Anz_Level=i-1;
      OldCell=n=Suche_Zelle(K,OldCell);
      Anz_Level=l;
      /**/
      /* Wenn nicht gefunden, dann entarteter Knoten */
      if(!n)
      {
        L.I=(Index_Typ *)realloc(L.I,(++L.Anz_Index)*sizeof(Index_Typ));
        L.I[L.Anz_Index-1]=j;
      }
    }
    /* Debug */
    printf("%d entartete Knoten\n",L.Anz_Index);
    /**/
    /* Nun haben wir eine vollständige Liste der entarteten Knoten für
       diesen Verfeinerungslevel. Nun suchen wir die Zellen, die 
       mindestens einen dieser Knoten enthalten.
    */
    /* Für alle Zellen des aktuellen Levels */
    for(j=1;j<=MG[i-1].Cells.Anz_Cells;j++)
    {
      n=0;
      /* Für jeden Knotenindex dieser Zelle */
      for(k=0;k<NpC;k++)
      /* Suche, ob in Liste */
      {
        for(l=0;l<L.Anz_Index;l++)
        {
          if (*(MG[i-1].Cells.Vert+j*NpC+k)==L.I[l])
          {
            n=1;
            break;
          }
        }
        if(n==1) break;
      }
      if(n==1)
      {
	CL[i-1].I=(Index_Typ *) realloc(CL[i-1].I,(++(CL[i-1].Anz_Index))
					*sizeof(Index_Typ));
	CL[i-1].I[CL[i-1].Anz_Index-1]=j;
      }
    }
    /* Debug */
    printf("%d entartete Zellen\n",CL[i-1].Anz_Index);
    /**/
  }
  free(L.I);
  /* Debug */
  printf("Startup_Suchen beendet!\n");
  /**/
}

void Cleanup_Suchen(void)
{
  Zaehler_Typ i;
  for(i=0;i<Anz_Level-1;i++) Cleanup_Cells(&(MG[i].Cells));
  Cleanup_Gitter(&(MG[Anz_Level-1]));
  Anz_Level=0;
  free(MG);
  for(i=0;i<Anz_Level;i++) 
  {
    CL[i].Anz_Index=0;
    if(CL[i].I!=NULL)
    {
      free(CL[i].I);
      CL[i].I=NULL;
    }
  }
}

Index_Typ Append_Tracers(Output_Typ *O,Anz_Typ Anz)
{
   Anz_Typ Neu,i;
   Neu=O->Tracers.Anz_Tracer+Anz;
   for(i=0;i<DIM;i++)
     O->Tracers.K[i]=(Daten_Typ *)realloc(O->Tracers.K[i],Neu*sizeof(Daten_Typ));
   O->Tracers.Farbe=(Daten_Typ *)realloc(O->Tracers.Farbe,Neu*sizeof(Daten_Typ));
   O->Tracers.Expire=(Daten_Typ *)realloc(O->Tracers.Expire,Neu*sizeof(Daten_Typ));
   O->Tracers.InCell=(Index_Typ *)realloc(O->Tracers.InCell,Neu*sizeof(Index_Typ));
   O->Tracers.B_Coords=(B_Coord_Typ *)realloc(O->Tracers.B_Coords,Neu*sizeof(B_Coord_Typ));
   for(i=0;i<Num_Field;i++)
     O->Field[i]=(Daten_Typ *)realloc(O->Field[i],Neu*sizeof(Daten_Typ));
   Neu=O->Tracers.Anz_Tracer;
   O->Tracers.Anz_Tracer+=Anz;
   return Neu;
}

void Clear_Tracer(Output_Typ *O,Index_Typ Pos)
{
  Index_Typ i,j;
  O->Tracers.Anz_Tracer--;
  for(i=Pos;i<O->Tracers.Anz_Tracer;i++)
  {
    for(j=0;j<DIM;j++)
      O->Tracers.K[j][i]=O->Tracers.K[j][i+1];
    O->Tracers.Farbe[i]=O->Tracers.Farbe[i+1];
    O->Tracers.Expire[i]=O->Tracers.Expire[i+1];
    O->Tracers.InCell[i]=O->Tracers.InCell[i+1];
    O->Tracers.B_Coords[i]=O->Tracers.B_Coords[i+1];
    for(j=0;j<Num_Field;j++)
      O->Field[j][i]=O->Field[j][i+1];
  }
}

#ifdef PART3D
/* muß noch verbessert werden */
void Calc_Barc(Gitter_Typ *G,Tracer_Typ *T,Index_Typ P)
{
  Anz_Typ i;
  for(i=0;i<NpC;i++)T->B_Coords[P].L[i]=1.0/(Daten_Typ)NpC;
}
#else

int Barc(Daten_Typ X1,Daten_Typ X2,Daten_Typ X3,Daten_Typ Y1,Daten_Typ Y2,
  Daten_Typ Y3,Daten_Typ X,Daten_Typ Y,Daten_Typ *L)
{
  Zaehler_Typ i;
  Daten_Typ Delta,A[3],B[3],C[3];
  Delta=(X3-X2)*(Y1-Y2)-(Y3-Y2)*(X1-X2);
  A[0]=(X2*Y3-X3*Y2)/Delta;
  B[0]=(Y2-Y3)/Delta;
  C[0]=(X3-X2)/Delta;
  A[1]=(X3*Y1-X1*Y3)/Delta;
  B[1]=(Y3-Y1)/Delta;
  C[1]=(X1-X3)/Delta;
  A[2]=(X1*Y2-X2*Y1)/Delta;
  B[2]=(Y1-Y2)/Delta;
  C[2]=(X2-X1)/Delta;
  for (i=0;i<3;i++) L[i]=A[i]+B[i]*X+C[i]*Y;
  for(i=0;i<3;i++)
  { 
    if(L[i]<=-Small_EPS) return 0;
  }
  return 1;
}

Daten_Typ Bounce(const Daten_Typ I)
{
  if(I<=0.0) return 0.0;
  if(I>=1.0) return 1.0;
  return I;
}

void Norm_Barc(Daten_Typ L[NpC])
{
  Zaehler_Typ i;
  Daten_Typ Erg;
  for(Erg=0.0,i=0;i<NpC;i++) Erg+=(L[i]=Bounce(L[i]));
  if(Erg<=Small_EPS) for(i=0;i<NpC;i++) L[i]=1.0/NpC;
  else for(i=0;i<NpC;i++) L[i]/=Erg;
}

void Calc_Barc(Gitter_Typ *G,Tracer_Typ *T,Index_Typ P)
{
  Index_Typ C;
  Anz_Typ NC;
  Zaehler_Typ i;
  Daten_Typ Alpha,Beta,Alph,Bet,X[4],Y[4],L[3];
  C=T->InCell[P];
  for(i=0;i<4;i++)
  {
    X[i]=G->Nodes.K[0][G->Cells.Vert[i+4*C]];
    Y[i]=G->Nodes.K[1][G->Cells.Vert[i+4*C]];
  }
  if (Is_Regular(C,G)==REG_REGULAR)
  {
    if(My_Abs(X[1]-X[0])>EPS) Alpha=(T->K[0][P]-X[0])/(X[1]-X[0]);
    else Alpha=(T->K[1][P]-Y[0])/(Y[1]-Y[0]);
    Alph=1.0-Alpha;
    if(My_Abs(X[0]-X[3])>EPS) Beta=(T->K[0][P]-X[0])/(X[3]-X[0]);
    else Beta=(T->K[1][P]-Y[0])/(Y[3]-Y[0]);
    Bet=1.0-Beta;
    T->B_Coords[P].L[0]=Alph*Bet;
    T->B_Coords[P].L[1]=Alpha*Bet;
    T->B_Coords[P].L[2]=Alpha*Beta;
    T->B_Coords[P].L[3]=Alph*Beta;  
  }
  else
  {
    NC=0;
    for(i=0;i<4;i++) T->B_Coords[P].L[i]=0.0;
    if(Barc(X[0],X[1],X[2],Y[0],Y[1],Y[2],T->K[0][P],T->K[1][P],&L[0]))
    {
      T->B_Coords[P].L[0]+=L[0];
      T->B_Coords[P].L[1]+=L[1];
      T->B_Coords[P].L[2]+=L[2];
      NC++;
    }
    if(Barc(X[0],X[2],X[3],Y[0],Y[2],Y[3],T->K[0][P],T->K[1][P],&L[0]))
    {
      T->B_Coords[P].L[0]+=L[0];
      T->B_Coords[P].L[2]+=L[1];
      T->B_Coords[P].L[3]+=L[2];
      NC++;
    }
    if(Barc(X[1],X[2],X[3],Y[1],Y[2],Y[3],T->K[0][P],T->K[1][P],&L[0]))
    {
      T->B_Coords[P].L[1]+=L[0];
      T->B_Coords[P].L[2]+=L[1];
      T->B_Coords[P].L[3]+=L[2];
      NC++;
    }
    if(Barc(X[0],X[1],X[3],Y[0],Y[1],Y[3],T->K[0][P],T->K[1][P],&L[0]))
    {
      T->B_Coords[P].L[0]+=L[0];
      T->B_Coords[P].L[1]+=L[1];
      T->B_Coords[P].L[3]+=L[2];
      NC++;
    }
    if (NC) for(i=0;i<4;i++) T->B_Coords[P].L[i]/=NC;
  }
  Norm_Barc(T->B_Coords[P].L);
}
#endif

void Gen_Tracers(Output_Typ *O,Block_Typ *B)
{
  Index_Typ Offset,OldCell;
  Zaehler_Typ i,j,l,Lev[DIM+1];
  Daten_Typ H[DIM],P[DIM];
  for(i=0,j=1;i<DIM;i++)
    j*=B->N[i];
  l=j;
  Offset=Append_Tracers(O,l);
  for(i=0;i<DIM;i++)
    H[i]=(B->P[1][i]-B->P[0][i])/B->N[i];
 
  Lev[0]=1;
  for(i=1;i<=DIM;i++) Lev[i]=Lev[i-1]*B->N[i-1];
  OldCell=0;
  for(i=0;i<l;i++)
    {
      for(j=0;j<DIM;j++)
	O->Tracers.K[j][Offset]=B->P[0][j]+(0.5+((i%Lev[j+1])/Lev[j]))*H[j];
      O->Tracers.Farbe[Offset]=B->Farbe;
      O->Tracers.Expire[Offset]=B->Haltbarkeit;
      for(j=0;j<DIM;j++)P[j]=O->Tracers.K[j][Offset];
      OldCell=O->Tracers.InCell[Offset]=Suche_Zelle(P,OldCell);
      if(!O->Tracers.InCell[Offset]) Clear_Tracer(O,Offset);
      else
	{
	  Calc_Barc(MG+Anz_Level-1,&(O->Tracers),Offset);
	  Offset++;
	}
    }
}

void Bewege_Tracer(Output_Typ *O,Input_Typ *I,Daten_Typ t)
{
  Zaehler_Typ i,j,k,pos,stop;
  Daten_Typ Dummy,P[DIM];
  stop=O->Tracers.Anz_Tracer;
  for(i=0,pos=0;i<stop;i++)
  {
    for(k=0;k<DIM;k++)
      {
	for(j=0,Dummy=0.0;j<NpC;j++)
	  Dummy+=O->Tracers.B_Coords[pos].L[j]*
	    I->V[k][I->Gitter.Cells.Vert[j+NpC*O->Tracers.InCell[pos]]];
	P[k]=O->Tracers.K[k][pos]+=t*Dummy;
      }
    O->Tracers.InCell[pos]=Suche_Zelle(P,O->Tracers.InCell[pos]);
    O->Tracers.Expire[pos]-=t;
    if((!O->Tracers.InCell[pos])||(O->Tracers.Expire[pos]<-EPS)) Clear_Tracer(O,pos);
    else
    {
      Calc_Barc(&(I->Gitter),&(O->Tracers),pos);
      pos++;
    }
  }
}

void Calculate_Values(Output_Typ *O,Input_Typ *I)
{
  Zaehler_Typ i,j,l;
  Daten_Typ Dummy;
  for(l=0;l<Num_Field;l++)
    {
      for(i=0;i<O->Tracers.Anz_Tracer;i++)
	{
	  for(j=0,Dummy=0.0;j<NpC;j++)
	    Dummy+=O->Tracers.B_Coords[i].L[j]*
	      I->Field[l][I->Gitter.Cells.Vert[j+NpC*O->Tracers.InCell[i]]];
	  O->Field[l][i]=Dummy;
	}
    }
}

const Daten_Typ DELTA=0.5E-2; 
void Schreibe_Output(const char *Name,Output_Typ *O,Gitter_Typ *F)
{
  FILE *f;
  Anz_Typ A,k;
 #ifdef MOVBC
  Anz_Typ i,j;
 #ifdef PART3D
 #define ANZ_POLY 6 
 const Index_Typ Index[3][2][ANZ_POLY*2]=
 {
   {
     {0,1,10,11,12,13,14,15,16,19,20,21},
     {2,3,4,5,6,7,8,9,17,18,22,23}
   },
   {
     {1,2,5,6,9,10,13,14,20,21,22,23},
     {0,3,4,7,8,11,12,15,16,17,18,19}
   },
   {
     {0,1,2,3,4,5,14,15,16,17,20,23},
     {6,7,8,9,10,11,12,13,18,19,21,22}
   }
 };
 #else
 #define ANZ_POLY 1 
 const Index_Typ Index[2][2][ANZ_POLY*2]=
 {
   {
     {0,3},
     {1,2}
   },
   {
     {0,1},
     {2,3}
   }
 };
 #endif
 Daten_Typ K[DIM],P[DIM][ANZ_POLY*4];
 #endif
  f=fopen(Name,"w");
  Schreibe_Gitter(f,&(O->Gitter));
#ifdef MOVBC
  /* Test (Materialien) */
  fprintf(f,"material 2 0\nmat1\nmat2\n");
  A=O->Gitter.Cells.Anz_Cells;
  for(i=0,j=1;i<A;i++) Schreibe_Index(f,1,&j);
  /* Test (Polygone) */
  fprintf(f,"polygons\n");
  A=F->Nodes.Anz_Nodes;
  for(i=1;i<=A;i++)
    {
    for(j=0;j<DIM;j++) K[j]=F->Nodes.K[j][i];
    if(Test_MovBC(K))
      {
	/* Ich baue zuerst eine Liste der Koordinaten
	   der einzelnen Polygone auf */
	for(k=0;k<DIM;k++)
	  {
	    for(j=0;j<ANZ_POLY*2;j++) 
	      {
	        P[k][Index[k][0][j]]=K[k]-DELTA;
	        P[k][Index[k][1][j]]=K[k]+DELTA;
	      } 
	  }
	/* Nun gebe ich diese Polygone aus */
	for(j=0;j<ANZ_POLY;j++)
	  {
	    fprintf(f,"2 4\n");
	    for(k=0;k<DIM;k++) Schreibe_Daten(f,4,P[k]+j*4);
	    if(DIM<3) Schreibe_Daten(f,4,NULL);
	  }
      }
    }
  fprintf(f,"endpoly\n");
  /**/
#endif
  fprintf(f,"tracers ");
  A=O->Tracers.Anz_Tracer;
  Schreibe_Anz(f,A);
  for(k=0;k<DIM;k++) Schreibe_Daten(f,A,O->Tracers.K[k]);
  if(DIM<3) Schreibe_Daten(f,A,NULL);
  for(k=0;k<Num_Field;k++)
  {
    fprintf(f,"%s\n",Tracer_Names[k]);
    Schreibe_Daten(f,A,O->Field[k]);
  }
  /* Das Feld "Farbe" faellt aus der Rolle, da es nicht aus den
     Simulationsdateien interpoliert wird und deshalb gesondert
     behandelt werden muß.
  */
  fprintf(f,"Farbe\n");
  Schreibe_Daten(f,A,O->Tracers.Farbe);
  fprintf(f,"endtrace\n");
  fprintf(f,"probtime ");
  Schreibe_Daten(f,1,&Global_Time);
  fprintf(f,"endgmv\n");
  fclose(f);
}

void Cleanup_Output(Output_Typ *O)
{
  Zaehler_Typ i;
  /* Gitter wird schon bei Cleanup_Suchen freigegeben! */
  for(i=0;i<DIM;i++) free(O->Tracers.K[i]);
  free(O->Tracers.Farbe);
  free(O->Tracers.Expire);
  free(O->Tracers.InCell);
  free(O->Tracers.B_Coords);
  O->Tracers.Anz_Tracer=0;
  for(i=0;i<Num_Field;i++) free(O->Field[i]);
}

Zaehler_Typ Modul_Input,Modul_Output;
Anz_Typ Anz_Zeitschritte,Anz_Block;
Daten_Typ Zeitschritt;
Block_Typ *Blocks;

void Lese_Parameter(const char *Name)
{
  FILE *f;
  Zaehler_Typ i,j;
  f=fopen(Name,"r");
  if(f==NULL) File_Error(Name);
  Lese_Anz(f,&Anz_Zeitschritte);Lese_Daten(f,1,&Zeitschritt);
  Lese_Anz(f,&Modul_Input);Lese_Anz(f,&Modul_Output);
  Lese_Anz(f,&Anz_Block);
  Blocks=(Block_Typ *)calloc(Anz_Block,sizeof(Block_Typ));
  for(i=0;i<Anz_Block;i++)
  {
    for(j=0;j<2;j++)
    {
      Lese_Daten(f,DIM,Blocks[i].P[j]);
    }
    Lese_Index(f,DIM,Blocks[i].N);
    Lese_Anz(f,&(Blocks[i].Modul));
    Lese_Anz(f,&(Blocks[i].Offset));
    Lese_Anz(f,&(Blocks[i].Start));
    Lese_Anz(f,&(Blocks[i].Stop));
    Lese_Daten(f,1,&(Blocks[i].Farbe));
    Lese_Daten(f,1,&(Blocks[i].Haltbarkeit));
  }
}

int main(int argc,char *argv[])
{
  Zaehler_Typ In_Pos,Out_Pos,T_Pos,i;
  char * Praefix;
  FILE *f;
  if (argc>1) Praefix=argv[1];
  else Praefix="t";
#ifdef PART3D
  printf("GMVPT-3D");
#else
  printf("GMVPT-2D");
#endif
  printf(" Version 1.1");
#ifdef MOVBC
  printf(" mit bewegtem Rand,");
#endif
#ifdef BOUSS
  puts(" angepasst fuer BOUSS...");
#else
  puts(" angepasst fuer PP2D,CP2D,...");
#endif
  Lese_Parameter("#data/gmvpt.dat");
  In_Pos=0;
  Out_Pos=0;
  for(i=0;i<DIM;i++) O.Tracers.K[i]=NULL;
  for(i=0;i<Num_Field;i++) O.Field[i]=NULL;
  O.Tracers.Anz_Tracer=0;
  O.Tracers.Farbe=NULL;
  O.Tracers.Expire=NULL;
  O.Tracers.InCell=NULL;
  O.Tracers.B_Coords=NULL;
  Startup_Suchen();
  f=fopen("#gmv/coarse.gmv","r");
  if(f==NULL) File_Error("#gmv/coarse.gmv");
  Lese_Gitter(f,&O.Gitter);
  fclose(f);
  for(T_Pos=0;T_Pos<Anz_Zeitschritte;T_Pos++)
  {
    if(T_Pos%Modul_Input==0)
    {
      do {
        In_Pos++;
        sprintf(NIn,"#gmv/u.%d.gmv",In_Pos);
      } while(!Datei_Existiert(NIn));
      if(T_Pos==0) Startup_Input(NIn,&I);
      else Update_Input(NIn,&I);
    }
    for(i=0;i<Anz_Block;i++)
    {
      if((T_Pos>=Blocks[i].Start)&&(T_Pos<=Blocks[i].Stop)&&
	 ((T_Pos+(Blocks[i].Modul-Blocks[i].Offset))%Blocks[i].Modul==0))
      {
        printf("Erzeuge Block %d ...\n",i+1); 
        Gen_Tracers(&O,Blocks+i);
      }
    }
    if(T_Pos%Modul_Output==0)
    {
      printf("Ausgabe der Tracer!\n");
      Out_Pos++;
      sprintf(NOut,"#gmv/%s.%d.gmv",Praefix,Out_Pos);
      Calculate_Values(&O,&I);
      Schreibe_Output(NOut,&O,&(I.Gitter));
    }
    printf("Bewege Tracer um einen Zeitschritt...\n");
    Global_Time+=Zeitschritt;
    Bewege_Tracer(&O,&I,Zeitschritt);

  }
  Cleanup_Input(&I);
  Cleanup_Output(&O);
  Cleanup_Suchen();
  return EXIT_SUCCESS;
}
