From 8fe367ecf2cbda79960a68a052fdb84916427c0e Mon Sep 17 00:00:00 2001
From: tradke <tradke@21a6bef8-4479-4f54-8f8d-0db94a2919ef>
Date: Wed, 13 Sep 2000 13:49:13 +0000
Subject: Importing latest stuff from development repository

git-svn-id: http://svn.cactuscode.org/arrangements/CactusExternal/FlexIO/trunk@2 21a6bef8-4479-4f54-8f8d-0db94a2919ef
---
 src/WildWriter.cc | 609 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 609 insertions(+)
 create mode 100644 src/WildWriter.cc

(limited to 'src/WildWriter.cc')

diff --git a/src/WildWriter.cc b/src/WildWriter.cc
new file mode 100644
index 0000000..b6bcf9c
--- /dev/null
+++ b/src/WildWriter.cc
@@ -0,0 +1,609 @@
+#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();
+}
-- 
cgit v1.2.3