#include <stdio.h>
#include <stdlib.h>
#include "WildWriter.hh"

//#define DISKBLOCKSIZE 8192

void WildWriterB::beginDataset(int numnodes,
			       int numcells,
			       IObase::DataType datacomponenttype,
			       int numdatacomponents,
			       int *componentsizes,
			       char *componentnames,
			       char componentnameseparator){
  int i;
  nnodes=numnodes;
  ncells=numcells;
  comptype=datacomponenttype;
  ncomp = numdatacomponents;
  compsizes.setSize(ncomp);
  if(componentsizes){
    for(i=0,totalcomp=1;i<ncomp;i++) {
      totalcomp+=compsizes[i]=1;
    }
  }
  else {
    for(i=0,totalcomp=1;i<ncomp;i++) {
      totalcomp+=compsizes[i]=1;
    }
  }
  compnames.setSize(strlen(componentnames)+2);
  compnames[0]=componentnameseparator;
  strcpy(compnames.getData()+1,componentnames);
  //  if(!buffer) buffer = new char[BLOCKSIZE];
}

void WildWriterB::beginNodes(){
  int dims[2]={3,0};
  dims[1]=nnodes;
  file.reserveChunk(IObase::Float32,2,dims);
  // set size of nodebuffer
  bufferindex=masterindex=0;
  segsize=BLOCKSIZE/(3*sizeof(float));
}

void WildWriterB::beginNodes(int numnodes){
  //setMode(WriteNodes);
  int dims[2]={3,0};
  nnodes=numnodes;
  dims[1]=nnodes;
  file.reserveChunk(IObase::Float32,2,dims);
  // set size of nodebuffer
  bufferindex=masterindex=0;
  segsize=BLOCKSIZE/(3*sizeof(float));
}
/*
  inline void addNode(float *coords){
    // attrib WildNode[%u] (should buffer internally)
    // but we can only chunk on Datasets :(
    // so collect data in 8k buffers
    checkMode(WriteNodes);
    float *copyptr=(float *)buffer;
    // copy into master buffer.
    int index=bufferindex*3;
    for(int i=0;i<3;i++) copyptr[index+i]=coords[i];
    bufferindex++;
    masterindex++;
    if(bufferindex>=segsize){ // dump the buffer to disk
      int origin[2]={0,0},dims[2]={3,0};
      dims[1]=segsize;
      origin[1]=masterindex-segsize;
      file.writeChunk(dims,origin,buffer);
      bufferindex=0; // reset the bufferindex for next seg
    }
  }
  */

void WildWriterB::endNodes(){
  // Flush buffers for final writeChunk()
  //setMode(Idle);
  int origin[2]={0,0},dims[2]={0,0};
  dims[0]=totalcomp;
  dims[1]=bufferindex;
  origin[1]=masterindex-bufferindex;
  // Flush buffers for final writeChunk()
  // dump the remainder of buffer
  file.writeChunk(dims,origin,buffer);
  bufferindex=0;
}

void WildWriterB::beginCells(){
  int dims[2]={9,0};
  dims[1]=ncells;
  file.reserveChunk(IObase::Int32,2,dims);
  bufferindex=masterindex=0;
  segsize=BLOCKSIZE/(9*sizeof(int));
}
void WildWriterB::beginCells(int numcells){
  int dims[2]={9,0};
  ncells=numcells;
  dims[1]=ncells;
  file.reserveChunk(IObase::Int32,2,dims);
  bufferindex=masterindex=0;
  segsize=BLOCKSIZE/(9*sizeof(int));
}
/*
  inline void addCell(int level,int *neighborlist){
    //checkMode(WriteCells);
    int *copyptr=(int *)buffer;
    int index=bufferindex*9;
    // copy into master buffer.
    copyptr[index+i]=level;
    for(int i=0;i<8;i++) 
      copyptr[index+i+1]=neighborlist[i];
    bufferindex++;
    masterindex++;
    if(bufferindex>=segsize){ // dump the buffer to disk
      int origin[2]={0,0},dims[2]={3,0};
      dims[1]=segsize;
      origin[1]=masterindex-segsize;
      file.writeChunk(dims,origin,buffer);
      bufferindex=0; // reset the bufferindex for next seg
    }
  }
  */

void WildWriterB::endCells(){
  // Flush buffers for final writeChunk()
  //setMode(Idle);
  int origin[2]={0,0},dims[2]={0,0};
  dims[0]=totalcomp;
  dims[1]=bufferindex;
  origin[1]=masterindex-bufferindex;
  // Flush buffers for final writeChunk()
  // dump the remainder of buffer
  file.writeChunk(dims,origin,buffer);
  bufferindex=0;
}
void WildWriterB::beginData(){
  int dims[2]={0,0};
  dims[0]=totalcomp;
  dims[1]=nnodes;
  file.reserveChunk(comptype,2,dims);
  bufferindex=masterindex=0;
  segsize=BLOCKSIZE/(totalcomp*IObase::sizeOf(comptype));
} 

void WildWriterB::beginData(IObase::DataType datacomponenttype,
			    int numdatacomponents,
			    int *componentsizes,
			    char *componentnames,
			    char componentnameseparator){
  int dims[2]={0,0};
  int i;
  comptype=datacomponenttype;
  ncomp = numdatacomponents;
  compsizes.setSize(ncomp);
  if(componentsizes){
    for(i=0,totalcomp=1;i<ncomp;i++)
      totalcomp+=compsizes[i]=componentsizes[i];
  }
  else {
    for(i=0,totalcomp=1;i<ncomp;i++)
      totalcomp+=compsizes[i]=1;
  }
  compnames.setSize(strlen(componentnames)+2);
  compnames[0]=componentnameseparator;
  strcpy(compnames.getData(),componentnames);
  dims[0]=totalcomp;
  dims[1]=nnodes;
  file.reserveChunk(comptype,2,dims);
  bufferindex=masterindex=0;
  segsize=BLOCKSIZE/(totalcomp*IObase::sizeOf(comptype));
}
/*
inline void WildWriterB::addData(void *data){
  //checkMode(WriteCells);
  if(comptype==IObase::Float32){
    float *copyptr=(float *)buffer;
    int index = bufferindex*totalcomp*sizeof(double);
    for(int i=0;i<totalcomp;i++) copyptr[index+i]=data[i];
  }
  else if(comptype==IObase::Float64){
    double *copyptr=(double *)buffer;
    int index = bufferindex*totalcomp*sizeof(double);
    for(int i=0;i<totalcomp;i++) copyptr[index+i]=data[i];
  }
  bufferindex++;
  masterindex++;
  if(bufferindex>=segsize){ // dump the buffer to disk
    int origin[2]={0,0},dims[2]={0,0};
    dims[0]=totalcomp;
    dims[1]=segsize;
    origin[1]=masterindex-segsize;
    file.writeChunk(dims,origin,buffer);
    bufferindex=0; // reset the bufferindex for next seg
  }
}
*/
void WildWriterB::endData(){
  int origin[2]={0,0},dims[2]={0,0};
  dims[0]=totalcomp;
  dims[1]=bufferindex;
  origin[1]=masterindex-bufferindex;
  // Flush buffers for final writeChunk()
  // dump the remainder of buffer
  file.writeChunk(dims,origin,buffer);
  bufferindex=0;
}

void WildWriterB::endDataset() {
  // flush all buffers (not really necessary in "dumb mode")
  // and make certain ncells & nnodes have been written
  // delete the master buffer if dynamically allocated
}

//#include "declare_extern.h"

void WriteNodes(char *filename){
  int i,j,k,d,cell_id,vertex_id;
  // FILE     *fp;
  // char     output_file[16]; 
  treenode *gp;
  // Later On, we'll need to keep the file open persistently
  IObase *iofile = new IEEEIO(filename,IObase::Write);
  WildWriterB *file = new WildWriterB(*iofile);

  for ( d=0 ; d<=total_depth ; d++ )
    for ( i=0 ; i<=(x[d]-1) ; i++ ) {
      gp = start_x[d][i];
      do {
	if ((gp->xyz[0] < xmax-e)&&
	    (gp->xyz[1] < ymax-e)&&
	    (gp->xyz[2] < zmax-e)){
	  cell_id++;
	  if (connectivity) {
	    if (gp->sib[1] != NULL)
	      if ((gp->sib[1]->xyz[0] < xmax-e)&&
		  (fabs(gp->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e))
		cell_id++;

	    if (gp->sib[3] != NULL)
	      if ((gp->sib[3]->xyz[1] < ymax-e)&&
		  (fabs(gp->sib[3]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e))
		cell_id++;

	    if ((gp->sib[1] != NULL)&&(gp->sib[3] != NULL))
	      if (gp->sib[3]->sib[1] != NULL)
		if ((gp->sib[3]->sib[1]->xyz[0] < xmax-e)&&
		    (gp->sib[3]->sib[1]->xyz[1] < ymax-e)&&
		    (fabs(gp->sib[3]->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e)&&
		    (fabs(gp->sib[3]->sib[1]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e))
		  cell_id++;

	    if (gp->sib[5] != NULL)
	      if ((gp->sib[5]->xyz[2] < zmax-e)&&
		  (fabs(gp->sib[5]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		cell_id++;

	    if (gp->sib[5] != NULL)
	      if (gp->sib[5]->sib[1] != NULL)
		if ((gp->sib[5]->xyz[2] < zmax-e)&&
		    (gp->sib[5]->sib[1]->xyz[0] < xmax-e)&&
		    (fabs(gp->sib[5]->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e)&&
		    (fabs(gp->sib[5]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		  cell_id++;

	    if (gp->sib[5] != NULL)
	      if (gp->sib[5]->sib[3] != NULL)
		if ((gp->sib[5]->xyz[2] < zmax-e)&&
		    (gp->sib[5]->sib[3]->xyz[1] < ymax-e)&&
		    (fabs(gp->sib[5]->sib[3]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e)&&
		    (fabs(gp->sib[5]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		  cell_id++;

	    if (gp->sib[5] != NULL)
	      if ((gp->sib[5]->sib[3] != NULL)&&(gp->sib[5]->sib[1] != NULL))
		if (gp->sib[5]->sib[3]->sib[1] != NULL)
		  if ((gp->sib[5]->sib[3]->sib[1]->xyz[0] < xmax-e)&&
		      (gp->sib[5]->sib[3]->sib[1]->xyz[1] < ymax-e)&&
		      (gp->sib[5]->sib[3]->sib[1]->xyz[2] < zmax-e)&&
		      (fabs(gp->sib[5]->sib[3]->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e)&&
		      (fabs(gp->sib[5]->sib[3]->sib[1]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e)&&
		      (fabs(gp->sib[5]->sib[3]->sib[1]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		    cell_id++;
	  }
	}
	gp = gp->sib[1];
      } while (gp != NULL);
    }
  for ( d=0 ; d<=total_depth ; d++ )
    for ( i=0 ; i<=(x[d]-1) ; i++ ) {
      gp = start_x[d][i];

      do {
	if ((gp->xyz[0] < xmax-e)&&
	    (gp->xyz[1] < ymax-e)&&
	    (gp->xyz[2] < zmax-e)) {
	  if (d == 0) {
	    gp->hex_id[0] = vertex_id+1;
	    gp->hex_id[1] = vertex_id+2;
	    gp->hex_id[2] = vertex_id+3;
	    gp->hex_id[3] = vertex_id+4;
	    gp->hex_id[4] = vertex_id+5;
	    gp->hex_id[5] = vertex_id+6;
	    gp->hex_id[6] = vertex_id+7;
	    gp->hex_id[7] = vertex_id+8;
	    vertex_id += 8;
	  }
	  else {
	    gp->hex_id[0] = gp->par->hex_id[gp->par_grid_point];
	    gp->hex_id[1] = vertex_id+1;
	    gp->hex_id[2] = vertex_id+2;
	    gp->hex_id[3] = vertex_id+3;
	    gp->hex_id[4] = vertex_id+4;
	    gp->hex_id[5] = vertex_id+5;
	    gp->hex_id[6] = vertex_id+6;
	    gp->hex_id[7] = vertex_id+7;
	    vertex_id += 7;
	  }
	}
	gp = gp->sib[1];
      } while (gp != NULL);
    }

  // Set the var names
  char str[256];
  str[0]='\0';
  for(i=0;i<variables;i++){
    char tmpstr[5];
    sprintf(tmpstr,"v%u,",i);
    strcat(str,tmpstr);
  }
  str[strlen(str)-1]='\0';
  /* Collected total sizes, now write UCD file header 
  fprintf(fp,"%d %d %d 0 0\n",vertex_id,cell_id,variables); */
  file->beginDataset(vertex_id, /* number of vertices */
		     cell_id, /* number of cells */
		     IObase::Float32, /* datatype for the data components */
		     variables, /* number of data components */
		     0,/* no array for componentsizes implies each is a scalar */
		     str,','); /* names */
  /* Now do the nodes */
  file->beginNodes();
  for ( d=0 ; d<=total_depth ; d++ )
    for ( i=0 ; i<=(x[d]-1) ; i++ ) {
      gp = start_x[d][i];
      do {
	if ((gp->xyz[0] < xmax-e)&&
	    (gp->xyz[1] < ymax-e)&&
	    (gp->xyz[2] < zmax-e)) {
	  float coords[3];
	  if (d == 0){
	    /* fprintf(fp,"%d %f %f %f\n",gp->hex_id[0], */
	    coords[0]=gp->xyz[0];
	    coords[1]=gp->xyz[1];
	    coords[2]=gp->xyz[2];
	    file->addNode(coords);
	  }
	  /* fprintf(fp,"%d %f %f %f\n",gp->hex_id[1], */
	  coords[0]=gp->xyz[0]+info[d][0];
	  coords[1]=gp->xyz[1];
	  coords[2]=gp->xyz[2];
	  file->addNode(coords);

	  /* fprintf(fp,"%d %f %f %f\n",gp->hex_id[2],*/
	  coords[0]=gp->xyz[0];
	  coords[1]=gp->xyz[1]+info[d][1];
	  coords[2]=gp->xyz[2];
	  file->addNode(coords);

	  /* fprintf(fp,"%d %f %f %f\n",gp->hex_id[3],*/
	  coords[0]=gp->xyz[0]+info[d][0];
	  coords[1]=gp->xyz[1]+info[d][1];
	  coords[2]=gp->xyz[2];
	  file->addNode(coords);

	  /*  fprintf(fp,"%d %f %f %f\n",gp->hex_id[4], */
	  coords[0]=gp->xyz[0];
	  coords[1]=gp->xyz[1];
	  coords[2]=gp->xyz[2]+info[d][2];
	  file->addNode(coords);

	  /* fprintf(fp,"%d %f %f %f\n",gp->hex_id[5], */
	  coords[0]=gp->xyz[0]+info[d][0];
	  coords[1]=gp->xyz[1];
	  coords[2]=gp->xyz[2]+info[d][2];
	  file->addNode(coords);

	  /* fprintf(fp,"%d %f %f %f\n",gp->hex_id[6], */
	  coords[0]=gp->xyz[0];
	  coords[1]=gp->xyz[1]+info[d][1];
	  coords[2]=gp->xyz[2]+info[d][2];
	  file->addNode(coords);

	  /*  fprintf(fp,"%d %f %f %f\n",gp->hex_id[7], */
	  coords[0]=gp->xyz[0]+info[d][0];
	  coords[1]=gp->xyz[1]+info[d][1];
	  coords[2]=gp->xyz[2]+info[d][2];
	  file->addNode(coords);
	}
	gp = gp->sib[1];
      } while (gp != NULL);
    }
  file->endNodes();
 
  file->beginCells(); 
  for ( cell_id=0,d=0 ; d<=total_depth ; d++ ){
    int cellvec[8],level;
    level=d;
    for ( i=0 ; i<=(x[d]-1) ; i++ ) {
      gp = start_x[d][i];
      do {
	if ((gp->xyz[0] < xmax-e)&&
	    (gp->xyz[1] < ymax-e)&&
	    (gp->xyz[2] < zmax-e))  {
	  cell_id++;

	  //fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
	  cellvec[0]=gp->hex_id[4];
	  cellvec[1]=gp->hex_id[5];
	  cellvec[2]=gp->hex_id[7];
	  cellvec[3]=gp->hex_id[6];
	  cellvec[4]=gp->hex_id[0];
	  cellvec[5]=gp->hex_id[1];
	  cellvec[6]=gp->hex_id[3];
	  cellvec[7]=gp->hex_id[2];
	  file->addCell(level,cellvec);
	  if (connectivity) {
	    if (gp->sib[1] != NULL)
	      if ((gp->sib[1]->xyz[0] < xmax-e)&&
		  (fabs(gp->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e))
		{
		  cell_id++;

		  // fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		  cellvec[0]=gp->hex_id[5];
		  cellvec[1]=gp->sib[1]->hex_id[4];
		  cellvec[2]=gp->sib[1]->hex_id[6];
		  cellvec[3]=gp->hex_id[7];
		  cellvec[4]=gp->hex_id[1];
		  cellvec[5]=gp->sib[1]->hex_id[0];
		  cellvec[6]=gp->sib[1]->hex_id[2];
		  cellvec[7]=gp->hex_id[3];
		  file->addCell(level,cellvec);
		}
	    
	    if (gp->sib[3] != NULL)
	      if ((gp->sib[3]->xyz[1] < ymax-e)&&
		  (fabs(gp->sib[3]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e))
		{
		  cell_id++;
		  //fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		  cellvec[0]= gp->hex_id[6];
		  cellvec[1]=gp->hex_id[7];
		  cellvec[2]=gp->sib[3]->hex_id[5];
		  cellvec[3]= gp->sib[3]->hex_id[4];
		  cellvec[4]=gp->hex_id[2];
		  cellvec[5]=gp->hex_id[3];
		  cellvec[6]= gp->sib[3]->hex_id[1];
		  cellvec[7]= gp->sib[3]->hex_id[0];
		  file->addCell(level,cellvec);
		}
	    
	    if ((gp->sib[1] != NULL)&&(gp->sib[3] != NULL))
	      if (gp->sib[3]->sib[1] != NULL)
		if ((gp->sib[3]->sib[1]->xyz[0] < xmax-e)&&
		    (gp->sib[3]->sib[1]->xyz[1] < ymax-e)&&
		    (fabs(gp->sib[3]->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e)&&
		    (fabs(gp->sib[3]->sib[1]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e))
		  {
		    cell_id++;
		    //fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		    cellvec[0]= gp->hex_id[7];
		    cellvec[1]= gp->sib[1]->hex_id[6];
		    cellvec[2]= gp->sib[3]->sib[1]->hex_id[4];
		    cellvec[3]= gp->sib[3]->hex_id[5];
		    cellvec[4]= gp->hex_id[3];
		    cellvec[5]= gp->sib[1]->hex_id[2];
		    cellvec[6]= gp->sib[3]->sib[1]->hex_id[0];
		    cellvec[7]= gp->sib[3]->hex_id[1];
		    file->addCell(level,cellvec);
		  }
	    
	    if (gp->sib[5] != NULL)
	      if ((gp->sib[5]->xyz[2] < zmax-e)&&
		  (fabs(gp->sib[5]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		{
		  cell_id++;
		  //fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		  cellvec[0]=gp->sib[5]->hex_id[0];
		  cellvec[1]=gp->sib[5]->hex_id[1];
		  cellvec[2]=gp->sib[5]->hex_id[3];
		  cellvec[3]=gp->sib[5]->hex_id[2];
		  cellvec[4]= gp->hex_id[4];
		  cellvec[5]=gp->hex_id[5];
		  cellvec[6]=gp->hex_id[7];
		  cellvec[7]= gp->hex_id[6];
		  file->addCell(level,cellvec);
		}
	    
	    if (gp->sib[5] != NULL)
	      if (gp->sib[5]->sib[1] != NULL)
		if ((gp->sib[5]->xyz[2] < zmax-e)&&
		    (gp->sib[5]->sib[1]->xyz[0] < xmax-e)&&
		    (fabs(gp->sib[5]->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e)&&
		    (fabs(gp->sib[5]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		  {
		    cell_id++;
		    //fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		    cellvec[0]=gp->sib[5]->hex_id[1];
		    cellvec[1]= gp->sib[5]->sib[1]->hex_id[0];
		    cellvec[2]= gp->sib[5]->sib[1]->hex_id[2];
		    cellvec[3]= gp->sib[5]->hex_id[3];
		    cellvec[4]=  gp->hex_id[5];
		    cellvec[5]= gp->sib[1]->hex_id[4];
		    cellvec[6]= gp->sib[1]->hex_id[6];
		    cellvec[7]=  gp->hex_id[7];
		    file->addCell(level,cellvec);
		  }

	    if (gp->sib[5] != NULL)
	      if (gp->sib[5]->sib[3] != NULL)
		if ((gp->sib[5]->xyz[2] < zmax-e)&&
		    (gp->sib[5]->sib[3]->xyz[1] < ymax-e)&&
		    (fabs(gp->sib[5]->sib[3]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e)&&
		    (fabs(gp->sib[5]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		  {
		    cell_id++;

		    //fprintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		    cellvec[0]=gp->sib[5]->hex_id[2];
		    cellvec[1]=gp->sib[5]->hex_id[3];
		    cellvec[2]=gp->sib[5]->sib[3]->hex_id[1];
		    cellvec[3]=gp->sib[5]->sib[3]->hex_id[0];
		    cellvec[4]=gp->hex_id[6];
		    cellvec[5]=gp->hex_id[7];
		    cellvec[6]=gp->sib[3]->hex_id[5];
		    cellvec[7]=gp->sib[3]->hex_id[4]; 
		    file->addCell(level,cellvec);
		  }
	    
	    if (gp->sib[5] != NULL)
	      if ((gp->sib[5]->sib[3] != NULL)&&(gp->sib[5]->sib[1] != NULL))
		if (gp->sib[5]->sib[3]->sib[1] != NULL)
		  if ((gp->sib[5]->sib[3]->sib[1]->xyz[0] < xmax-e)&&
		      (gp->sib[5]->sib[3]->sib[1]->xyz[1] < ymax-e)&&
		      (gp->sib[5]->sib[3]->sib[1]->xyz[2] < zmax-e)&&
		      (fabs(gp->sib[5]->sib[3]->sib[1]->xyz[0]-gp->xyz[0]) <= info[d][0]+info[d][0]+e)&&
		      (fabs(gp->sib[5]->sib[3]->sib[1]->xyz[1]-gp->xyz[1]) <= info[d][1]+info[d][1]+e)&&
		      (fabs(gp->sib[5]->sib[3]->sib[1]->xyz[2]-gp->xyz[2]) <= info[d][2]+info[d][2]+e))
		    {
		      cell_id++;
		      //rintf(fp,"%d 1 hex %d %d %d %d %d %d %d %d\n",cell_id,
		      cellvec[0]=gp->sib[5]->hex_id[3];
		      cellvec[1]=gp->sib[5]->sib[1]->hex_id[2];
		      cellvec[2]=gp->sib[5]->sib[3]->sib[1]->hex_id[0];
		      cellvec[3]=gp->sib[5]->sib[3]->hex_id[1];
		      cellvec[4]=gp->hex_id[7];
		      cellvec[5]=gp->sib[1]->hex_id[6];
		      cellvec[6]=gp->sib[3]->sib[1]->hex_id[4];
		      cellvec[7]=gp->sib[3]->hex_id[5];
		      file->addCell(level,cellvec);
		    }
	  }
	}
	gp = gp->sib[1];
      } while (gp != NULL);
    }
  }
  file->endCells(); // Done with writing Cells
  //rintf(fp,"%d",variables);
  //for ( i=1 ; i<=(variables-1) ; i++ )
  //  fprintf(fp," 1");
  // fprintf(fp," 1\n");
  file->beginData(); // Start writing Data
  //for ( i=1 ; i<=variables ; i++ )
  //  fprintf(fp,"v%d, 1b/in**2\n",i);
  for ( d=0 ; d<=total_depth ; d++ )
    for ( i=0 ; i<=(x[d]-1) ; i++ ) {
      gp = start_x[d][i];
      do {
	if ((gp->xyz[0] < xmax-e)&&
	    (gp->xyz[1] < ymax-e)&&
	    (gp->xyz[2] < zmax-e)) {
	  if (d == 0)  {  // The toplevel 
	    for ( k=0 ; k<=7 ; k++ ) {
	      /* fprintf(fp,"%d ",gp->hex_id[k]);
	      for ( j=0 ; j<=(variables-2) ; j++ )
			fprintf(fp,"%f ",gp->u[k][t[d][0]][j]);
	      fprintf(fp,"%f\n",gp->u[k][t[d][0]][variables-1]);*/
	      file->addData(gp->u[k][t[d][0]]);
	    }
	  }
	  else {
	    for ( k=1 ; k<=7 ; k++ ) {
	      /*fprintf(fp,"%d ",gp->hex_id[k]);
	      for ( j=0 ; j<=(variables-2) ; j++ )
			fprintf(fp,"%f ",gp->u[k][t[d][0]][j]);
	      fprintf(fp,"%f\n",gp->u[k][t[d][0]][variables-1]);*/
	      file->addData(gp->u[k][t[d][0]]);
	    }
	  }
	}
	gp = gp->sib[1];
      } while (gp != NULL);
    }
  file->endData();
  file->endDataset();
}