sandbox/sander/output_htg.h

    #include "output_pvd.h"
    #include "utils.h"

    This field enables the export of AMR-meshes into the HyperTreeGrid Dataformat. It works on Quad-/Octrees with and without MPI. The HTG file format is still under development so compatibility with future Paraview versions might be limited. Advantage of the HTG file format is the ~7x smaller size due to implicit point location. The data is stored under the path variable, a time series can be opend using the corresponding collection file (.pvd) with name “output_.pvd”

    output_htg() is a wrapper function that additionally creates a collection file pointing to different time steps of an htg simulation. If MPI is used the HTG is written with collective MPI_File_IO. If no MPI is used standard POSIX fopen/fclose operations are used.

    Be sure to used the updated, MPI_IO compatible output_pvd.h!

    Example Usage

    event report(i+=10){
      char path[]="htg"; // no slash at the end!!
      char prefix[80];
      sprintf(prefix, "data_%06d", i);
      output_htg((scalar ) {cs, d,p},(vector ){u}, path, prefix, i, t);
    }

    The Collection File Name includes the used path. For the example above the Collection File is called “output_htg”. This allows to run a parameter study based on the same binary if you include the current case configuration in the path name. (e.g. Level, Reynolds number, velocity, acceleration, …)

    Flags

    #define HEADER_MIN_MAX_VAL 1

    If this flag is set the corresponding min/max values for each values are included in the header. This is probably not needed, the adavantage is unknown. It probably increases the posprocessing speed. Defaults to 1 (true).

    #define WRITE_HTG_LEVEL 1

    If this flag is set the level of each cell is exported as well and selectable during postprocessing. Defaults to 0 (false).

    #define HTG_SPEED_STATS 1

    If this flag is set information about the write speed are posted (to stdout). Defaults to 0 (false).

    #define VTK_FILE_VERSION 20

    If this flag is set to 20, the created HyperTreeGrid has the VTKFile version 2.0. It is supported starting with Paraview Version 5.10. Defaults to 10 (VTKFile version 1.0)

    Known HTG Problems

    • HyperTreeGridToDualGrid stops working when advancing time step (Paraview related)
    • Contour Filter does not work on HTG with only one tree (like this exporter creates) (Paraview related)
    • x-z Axis are swapped (3D), x-y Axis are swapped (2D)

    Details of the parallel implementation

    The mapping of the Basilisk internal tree structure to the HyperTreeGrid structure is not trivial. The HTG structure is defined level wise. Multiple processors can write a single tree, if each processes writes its cells level wise after those from the preceeding process. The offsets and cells per level have to recalculated every output step. These offsets and cells per level are used to create a derived datatype (MPI_Type_indexed). Each process copys the values of all its nodes into a continous array. The MPI_Type_indexed datatype allows to write this continous array as a correclty spaced array to the file, which allows other processes to write their data in the spaces left. The HyperTreeGrid further needs the size of the following data in bytes in addition to the actual data.

     _| size_in_bytes_of_the_following_blob |_actual_data_scalar1 | size...blob |_ac...data |...

    It is therefore convenient to define a struct consisting of the an integer with the blob size and the array of the data to write. This can be done using an MPI_Type_struct datatype. While the struct contains the interger and a pointer to a continous space in memory the data on the file will be the integer followed by the correctly spaced MPI_Type_indexed.

    This enables using MPI_File_set_view and MPI_File_write_all (collective operation) for good IO-performance.

    The header and the tail is written only by process 0.

    The actual tree structure of the HTG ist defined using a bitmask. If a node has children its value is 1, else it is 0. This is done level wise as well. On the next level the value of each (still existing) node is described using 1 (is refined/has childen) or 0 (not refined/ does not have childen) as well. This is repeated for each level except the last level because these values are all 0 due to being the last level and are neglected.

    The problem arises when storing this bitmask, as it is only possible the write a byte (8 bits). If a process on a level does not have a number of nodes that is devisible by 8 (a whole byte), logical gaps in the bitmask appear which are not allowed. This problem arises immediately at level 0, as process 0 has one 1 bit (level zero has children), but writing a 1000 0000 byte is wrong as it already describes the values of the 7 following nodes as well, whose values are unknown to process 0. Therefore the describtor bits (stored as bytes (u_int8_t) in memory) have to be moved to the next process (or level) until a number divisible by 8 is resident on the process. A ‘wave’ of excess descriptor bytes is swept through all processes and levels til the last process of the last level who fills the missing bytes with 0.

    To allow this moving of upto 7 bytes a special array structure is chosen:

    |8_empty_byt|descr_byt_lvl0|7_empty_byt|descr_byt_lvl1|…|7_empty_byt| 0 8

    Each process has to space to accommodate upto 7 bytes from the preceeding process/level infront of the own descriptor bytes. Each process keeps track of how many bytes to send/ to recv/ and stay on a level and on which index the level begins taking the recvieved bytes into account. After this every process has a number of descriptor bytes that can be transformed into descriptor bits. This is actually done in the same array starting a index 0 which is guaranteed to not contain any data. In byte 0 the information of the next 8 relevant bytes is written as bits.

    With information about offset and length again a MPI_Type_indexed and a MPI_Type struct can be defined. This again is used to write the bitmask in a collective operation.

    Detail of the serial implementation

    Each process caches the data on each level in a contingues array which is written to disk sequentially with only one fwrite() funtion call. The Bitmask is written in “appened” mode as well resulting in an uncluttered file

    //~ #define HEADER_MIN_MAX_VAL 0
    //~ #define WRITE_HTG_LEVEL 1
    //~ #define HTG_SPEED_STATS 1
    //~ #define VTK_FILE_VERSION 20 // 2.0 
    
    #ifndef HEADER_MIN_MAX_VAL
    #define HEADER_MIN_MAX_VAL 1
    #endif
    #ifndef WRITE_HTG_LEVEL
    #define WRITE_HTG_LEVEL 0
    #endif
    #ifndef HTG_SPEED_STATS
    #define HTG_SPEED_STATS 0
    #endif
    #ifndef VTK_FILE_VERSION
    #define VTK_FILE_VERSION 10 // 1.0
    #endif
    
    void output_htg(scalar * list, vector * vlist, const char* path, char* prefix, int i, double t);
    
    #if _MPI
    void output_htg_data_mpiio(scalar * list, vector * vlist, MPI_File fp);
    #else
    void output_htg_data(scalar * list, vector * vlist, FILE * fp);
    #endif
     
    #if _MPI
    void output_htg(scalar * list, vector * vlist, const char* path, char* prefix, int i, double t) {
    	MPI_File fp;
    	
    	int ec;
    	char htg_name[80];  	  
    	sprintf(htg_name, "%s/%s.htg", path, prefix);  
    
    	ec = MPI_File_open(MPI_COMM_WORLD, htg_name, \
    		MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &fp);
      // Overwrite File 
    	if (ec == MPI_ERR_FILE_EXISTS) {
        printf("ERR, htg_name exists!\n");
    		// MPI_File_delete(htg_name, MPI_INFO_NULL);
    		MPI_File_open(MPI_COMM_WORLD, htg_name, \
    		  MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &fp);
        MPI_File_set_size(fp, 0);
    	}
    	
    	if(ec != MPI_SUCCESS){
    		printf("output_htg.h : %s could not be opened\n Does the Folder exist?\n", htg_name);
    		MPI_Abort(MPI_COMM_WORLD, 2);
    	}  
    	
    	output_htg_data_mpiio((scalar *) list,(vector *)vlist,fp);
    
    	MPI_File_close(&fp);
    	
      if(pid()==0){
    		bool firstTimeWritten = false;
    		char pvd_name[80];	  
    		sprintf(pvd_name,"output_%s.pvd",path);
    	
    		ec = MPI_File_open(MPI_COMM_SELF, pvd_name, \
    			MPI_MODE_RDWR, MPI_INFO_NULL, &fp);
          
    		if( (i == 0) ||  (ec == MPI_ERR_NO_SUCH_FILE) ) {
    			// Overwrite File 	
    			// Best implementation??		
    			MPI_File_open(MPI_COMM_SELF, pvd_name, \
    			  MPI_MODE_WRONLY|MPI_MODE_CREATE, MPI_INFO_NULL, &fp);			
    			MPI_File_set_size(fp, 0);
    			firstTimeWritten = true;
    		}
    		output_pvd_mpiio(htg_name, t, fp, firstTimeWritten);
    		
    		MPI_File_close(&fp);
    	}	
    	MPI_Barrier(MPI_COMM_WORLD);    
    }
    
    #else // no MPI
    
    void output_htg(scalar * list, vector * vlist, const char* path, char* prefix, int i, double t) {
    	FILE * fp ;
    
    	char htg_name[80];  	  
    	sprintf(htg_name, "%s/%s.htg", path, prefix);  
    
    	fp = fopen(htg_name, "w");	  
    	if(!fp){
    		printf("output_htg.h : %s could not be opened\n Does the Folder exist?\n", htg_name);
    		exit(1);
    	}
      
      output_htg_data((scalar *) list,(vector *) vlist, fp);
    
      fclose(fp);
    
      bool firstTimeWritten = false;
      char pvd_name[80];	  
      sprintf(pvd_name,"output_%s.pvd",path);
      fp = fopen(pvd_name, "r+");
      if( (i == 0) ||  (fp == NULL) ) {
        fp = fopen(pvd_name,"w");
        firstTimeWritten = true;
      }
      output_pvd(htg_name, t, fp, firstTimeWritten);
      fclose(fp);	
    }
    
    #endif
    
    
    #if _MPI
    
    #define Write2File(x) do { \
      (x); \
      MPI_File_write(fp,&buffer, strlen(buffer), MPI_CHAR, MPI_STATUS_IGNORE); \
      } while(0)
      
    void output_htg_data_mpiio(scalar * list, vector * vlist, MPI_File fp)
    {
    	#if defined(_OPENMP)
    		int num_omp = omp_get_max_threads();
    		omp_set_num_threads(1);
    	#endif
        
    	unsigned int vertices_local = 0;
    	unsigned int descBits_local = 0;  
      unsigned int vertices_local_pL[grid->maxdepth+1];  // pL = per Level 
      
      unsigned int descBits;
      unsigned int vertices;  
      unsigned int vertices_pL[grid->maxdepth+1];  
    	    
    	for (int lvl = 0; lvl < grid->maxdepth+1; ++lvl){
        vertices_local_pL[lvl] = 0;
    		foreach_level(lvl,serial,noauto) 
    			if(is_local(cell)) 
    				vertices_local_pL[lvl]++;
        
        vertices_local += vertices_local_pL[lvl];         
      }
      
      descBits_local = vertices_local - vertices_local_pL[grid->maxdepth]; // 
      
      MPI_Reduce(&vertices_local_pL[0], &vertices_pL[0], grid->maxdepth + 1,\
                  MPI_UNSIGNED, MPI_SUM, 0, MPI_COMM_WORLD);  
     
    #if HTG_SPEED_STATS
    	MPI_Barrier(MPI_COMM_WORLD);
    	double start = MPI_Wtime();
    #endif	
      MPI_Offset offset = 0;
    
    	int vertices_global_offset[grid->maxdepth+1];
      vertices_global_offset[0] = 0;
      unsigned int carryover = 0;
      for (int lvl = 0; lvl <= grid->maxdepth; ++lvl){
        MPI_Exscan (&vertices_local_pL[lvl], \
                    &vertices_global_offset[lvl], \
                    1, \
                    MPI_UNSIGNED, \
                    MPI_SUM, \
                    MPI_COMM_WORLD);    
        // Pass Offset to process 0 for Exscan on next level
        if (pid() == (npe() - 1)) {
          unsigned int next_offset;
          next_offset = vertices_global_offset[lvl] + vertices_local_pL[lvl];
          MPI_Ssend(&next_offset, 1, MPI_UNSIGNED, 0, 0, MPI_COMM_WORLD);            
        }
        if (pid() == 0) {
          vertices_local_pL[lvl] -= carryover; // temporarily remove offset added in previous level
          
          MPI_Recv(&carryover, 1, MPI_UNSIGNED, npe() - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);      
          
          if (lvl < grid->maxdepth) {
            vertices_local_pL[lvl+1] += carryover; // temporarily use array to carry over offset to next level
            vertices_global_offset[lvl+1] = carryover;
          }      
          else
            vertices = carryover;
          if (lvl == grid->maxdepth - 1)
            descBits = carryover;
        }
      }
      
      MPI_Bcast(&vertices, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);
      MPI_Bcast(&descBits, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);
    • Calculate min,max for xml header (prob. not needed?)
    #if HTG_SPEED_STATS     
      MPI_Barrier(MPI_COMM_WORLD);
      double t_stats_start = MPI_Wtime();
    #endif  
      double min_val[list_len(list)];
      double max_val[list_len(list)];
      {
        int i = 0;
        for (scalar s in list) {
    #if HEADER_MIN_MAX_VAL
          stats stat = statsf(s);
          min_val[i]=stat.min;
          max_val[i]=stat.max;
    #else
          min_val[i]=0.;
          max_val[i]=0.;
    #endif
          i++;
        }
    	}
      double min_val_v[vectors_len(vlist)];
      double max_val_v[vectors_len(vlist)];
      { 
        int i = 0;	
        for (vector v in vlist) {
    #if HEADER_MIN_MAX_VAL
          min_val_v[i]= HUGE;
          max_val_v[i]= -HUGE;
          foreach_dimension(){
            stats stat = statsf(v.x);	
            min_val_v[i] = min(stat.min,min_val_v[i]);
            max_val_v[i]=  max(stat.max,max_val_v[i]);			
          }
    #else
          min_val_v[i]= 0.;
          max_val_v[i]= 0.;
    #endif
          i++;
        }
    	}
    #if HTG_SPEED_STATS    
      MPI_Barrier(MPI_COMM_WORLD);
      double t_stats = MPI_Wtime() - t_stats_start ;
    	double t_header_start = MPI_Wtime();
    #endif  
    	char buffer[256];

    File Header vtkHyperTreeGrid

    	if (pid() == 0){
        //~ WriteHeader(&offset, infoStruct);
        int maj_v = VTK_FILE_VERSION/10, min_v = VTK_FILE_VERSION%10;
        
        Write2File(sprintf(buffer,\
          "<VTKFile %s version=\"%i.%i\" %s %s>\n",\
            "type=\"HyperTreeGrid\"", \
            maj_v, min_v , \
            "byte_order=\"LittleEndian\" ",\
            "header_type=\"UInt32\"")); 
      
      
      #if dimension==2
    		Write2File(sprintf(buffer,\
          "\t<HyperTreeGrid BranchFactor=\"2\" TransposedRootIndexing=\"0\" Dimensions=\"%d %d %d\">\n",\
          2,2,1));    
      #elif dimension==3
    		Write2File(sprintf(buffer,\
          "\t<HyperTreeGrid BranchFactor=\"2\" TransposedRootIndexing=\"0\" Dimensions=\"%d %d %d\">\n", 2,2,2));
      #endif
    		Write2File(sprintf(buffer,"\t\t<Grid>\n"));
      #if dimension==2
    		Write2File(sprintf(buffer, "\t\t\t<DataArray type=\"Float64\" Name=\"XCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",Y0, Y0+L0));
        Write2File(sprintf(buffer, "\t\t\t\t%g %g",Y0, Y0+L0));
        Write2File(sprintf(buffer,"\n\t\t\t</DataArray>\n"));
        Write2File(sprintf(buffer, "\t\t\t<DataArray type=\"Float64\" Name=\"YCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",X0, X0+L0));
        Write2File(sprintf(buffer, "\t\t\t\t%g %g",X0, X0+L0));
        Write2File(sprintf(buffer,"\n\t\t\t</DataArray>\n"));	  
    		Write2File(sprintf(buffer, "\t\t\t<DataArray type=\"Float64\" Name=\"ZCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",Z0, Z0+L0));		
    		Write2File(sprintf(buffer, "\t\t\t\t%g %g", 0., 0.));	  
      #elif dimension==3
    		Write2File(sprintf(buffer, "\t\t\t<DataArray type=\"Float64\" Name=\"XCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",\
            Z0, Z0+L0));
    		Write2File(sprintf(buffer, "\t\t\t\t%g %g",Z0, Z0+L0));
    		Write2File(sprintf(buffer,"\n\t\t\t</DataArray>\n"));
    		Write2File(sprintf(buffer, "\t\t\t<DataArray type=\"Float64\" Name=\"YCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",\
            Y0, Y0+L0));
    		Write2File(sprintf(buffer, "\t\t\t\t%g %g",Y0, Y0+L0));
    		Write2File(sprintf(buffer,"\n\t\t\t</DataArray>\n"));	  
    		Write2File(sprintf(buffer, "\t\t\t<DataArray type=\"Float64\" Name=\"ZCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",\
            X0, X0+L0));		
    		Write2File(sprintf(buffer, "\t\t\t\t%g %g",X0, X0+L0));	  
      #endif
    		Write2File(sprintf(buffer,"\n\t\t\t</DataArray>\n"));	  
    		Write2File(sprintf(buffer,"\t\t</Grid>\n"));
    		Write2File(sprintf(buffer,"\t\t<Trees>\n"));
        
    		// Tree Begin
        unsigned int byte_offset = 0;
      #if VTK_FILE_VERSION == 10
    		Write2File(sprintf(buffer,"\t\t\t<Tree Index=\"0\" NumberOfLevels=\"%d\" NumberOfVertices=\"%u\">\n",\
            grid->maxdepth + 1, vertices));   	
    		
        // Array Descriptor Bits    
        Write2File(sprintf(buffer,"\t\t\t\t<DataArray type=\"Bit\" Name=\"Descriptor\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"0\" RangeMax=\"1\" offset=\"%u\"/>\n",\
          descBits, byte_offset));        
        byte_offset += (descBits/8+1) * sizeof(u_int8_t) + sizeof(u_int32_t);
        
    		// Array NbVerticesByLevel
    		Write2File(sprintf(buffer,"\t\t\t\t<DataArray type=\"Int64\" Name=\"NbVerticesByLevel\" NumberOfTuples=\"%d\" format=\"ascii\" RangeMin=\"1\" RangeMax=\"%u\" >\n\t\t\t\t\t",\
            grid->maxdepth + 1, vertices_pL[grid->maxdepth]));
        		
        for (int lvl = 0; lvl <= grid->maxdepth; lvl++) {
    			Write2File(sprintf(buffer,"%u ", vertices_pL[lvl]));
    		}		
    		Write2File(sprintf(buffer,"\n\t\t\t\t</DataArray>\n"));
      #endif     
      
      #if VTK_FILE_VERSION == 20
        // Array Descriptor Bits    
        Write2File(sprintf(buffer,"\t\t\t\t<DataArray type=\"Bit\" Name=\"Descriptors\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"0\" RangeMax=\"1\" offset=\"%u\"/>\n",\
          descBits, byte_offset));        
        byte_offset += (descBits/8+1) * sizeof(u_int8_t) + sizeof(u_int32_t);
        
        Write2File(sprintf(buffer,"\t\t\t\t<DataArray type=\"Int64\" Name=\"NumberOfVerticesPerDepth\" NumberOfTuples=\"%d\" format=\"ascii\" RangeMin=\"1\" RangeMax=\"%u\" >\n\t\t\t\t\t",\
            grid->maxdepth + 1, vertices_pL[grid->maxdepth]));
              
        for (int lvl = 0; lvl <= grid->maxdepth; lvl++) {
          Write2File(sprintf(buffer,"%u ", vertices_pL[lvl]));
        }		
        Write2File(sprintf(buffer,"\n\t\t\t\t</DataArray>\n"));
        
        Write2File(sprintf(buffer,"\t\t\t\t<DataArray type=\"Int64\" Name=\"TreeIds\" NumberOfTuples=\"1\" format=\"ascii\" RangeMin=\"0\" RangeMax=\"0\" >\n"));
        Write2File(sprintf(buffer,"\t\t\t\t\t%i\n",0));
        Write2File(sprintf(buffer,"\t\t\t\t</DataArray>\n"));
        
        Write2File(sprintf(buffer,"\t\t\t\t<DataArray type=\"UInt32\" Name=\"DepthPerTree\" NumberOfTuples=\"1\" format=\"ascii\" RangeMin=\"%i\" RangeMax=\"%i\" >\n",\
        grid->maxdepth+1, grid->maxdepth+1));
        Write2File(sprintf(buffer,"\t\t\t\t\t%i\n",grid->maxdepth+1));
        Write2File(sprintf(buffer,"\t\t\t\t</DataArray>\n"));
                
        Write2File(sprintf(buffer,"\t\t</Trees>\n"));
      #endif  
    		// Cell Data Begin
    		Write2File(sprintf(buffer,"\t\t\t\t<CellData>\n"));
        
    #if WRITE_HTG_LEVEL    
    		Write2File(sprintf(buffer,"\t\t\t\t\t<DataArray type=\"UInt8\" Name=\"Level\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"0\" RangeMax=\"%i\" offset=\"%u\"/>\n",\
            vertices, grid->maxdepth, byte_offset));		
    		byte_offset += vertices * sizeof(u_int8_t) + sizeof(u_int32_t);
    #endif    
    		{
          int i = 0;
          for (scalar s in list)
          {			
            Write2File(sprintf(buffer,"\t\t\t\t\t<DataArray type=\"Float32\" Name=\"%s\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"%g\" RangeMax=\"%g\" offset=\"%u\"/>\n",\
                s.name, vertices, min_val[i], max_val[i], byte_offset));			    
            byte_offset += vertices * sizeof(float_t) + sizeof(u_int32_t);
            i++;
          }
        }
        {
          int i = 0;
          for (vector v in vlist)
          {
            char *vname = strtok(v.x.name, ".");
            Write2File(sprintf(buffer, "\t\t\t\t\t<DataArray type=\"Float32\" NumberOfComponents=\"%i\" Name=\"%s\" NumberOfTuples=\"%u\" format=\"appended\"  RangeMin=\"%g\" RangeMax=\"%g\" offset=\"%u\"/>\n",\
                3,vname, vertices, min_val_v[i],max_val_v[i],byte_offset));
            byte_offset += vertices * 3 * sizeof(float_t) + sizeof(u_int32_t);
            i++;
          } 
        }
    		// Cell Data End  		
    		Write2File(sprintf(buffer,"\t\t\t\t</CellData>\n"));
    #if VTK_FILE_VERSION == 10
        Write2File(sprintf(buffer,"\t\t\t</Tree>\n\t\t</Trees>\n"));
    #endif
    
        Write2File(sprintf(buffer,"\t</HyperTreeGrid>\n\t<AppendedData encoding=\"raw\">\n_"));
        MPI_Offset offset_tmp;
        MPI_File_get_position(fp, &offset_tmp);
        offset += offset_tmp;    
    	} // end pid() == 0
    	
    	MPI_Bcast(&offset, 1, MPI_OFFSET, 0, MPI_COMM_WORLD);
    #if HTG_SPEED_STATS  
      	double t_header = MPI_Wtime() - t_header_start;
    #endif

    Write Bitmask

      int cell_size;
      
      { 
        cell_size=sizeof(u_int8_t);
        int vertices_local_pL_offset[grid->maxdepth+1];

    Create an array large enough to hold data + spacing between the levels

        long length_w_spacing = descBits_local + 7*(grid->maxdepth+1) + 8;    
        u_int8_t* mask = (u_int8_t*)calloc( length_w_spacing, cell_size); //calloc needed?
        
        long index = 8; // start a 8
        for(int lvl=0; lvl < grid->maxdepth;++lvl) { // iterate through array and place lokal mask.
          vertices_local_pL_offset[lvl] = index;
          foreach_level(lvl,serial,noauto) {
            if (is_local(cell)) {
              mask[index++] = (u_int8_t)(!is_leaf(cell));
            }
          }
          index += 7; // make space for byte movement
        }
        vertices_local_pL_offset[grid->maxdepth] = index;
        
        assert(length_w_spacing > index);
      
        int vertices_local_pL_corr[grid->maxdepth+1];
        
        for (int lvl = 0; lvl < grid->maxdepth; ++lvl) {
          vertices_local_pL_corr[lvl] = vertices_local_pL[lvl];
        }
        
        
        int vertices_local_pL_offset_corr[grid->maxdepth]; // new offset
        long vertices_local_corr = 0;
        
        int num_from_prev = 0;
        u_int8_t tmp[7];
        for (int lvl = 0; lvl < grid->maxdepth; ++lvl) {
          for (int pe=0; pe < npe();  ++pe) {        
            int send_rank = pe;
            int recv_rank = (pe+1) % npe(); // next higher rank + loop
                     
            if (pid() == send_rank) {
              // work on prev send data          
              vertices_local_pL_corr[lvl] += num_from_prev;          
              int trg_position = vertices_local_pL_offset[lvl]-num_from_prev;
              vertices_local_pL_offset_corr[lvl] = trg_position;            
              for (int tmp_cnt = 0; tmp_cnt < num_from_prev; ++tmp_cnt)
              {
                // printf("pid: %i, lvl: %i trg_pos: %i, tmp_cnt: %i, lws: %li\n", pid(), lvl, trg_position, tmp_cnt,length_w_spacing); fflush(stdout);
                mask[trg_position + tmp_cnt] = tmp[tmp_cnt];
              }
              int num_to_next = (vertices_local_pL_corr[lvl]%8);
             
              vertices_local_pL_corr[lvl] -= num_to_next;
              
              // if last level and last process, append space 
              if ((lvl ==  grid->maxdepth -1 ) && (pid() == npe() -1) ) {
                vertices_local_pL_corr[lvl] += 8;            
              }
              
              vertices_local_corr += vertices_local_pL_corr[lvl];               
              
              // MPI_Ssend(&num_to_next,  1, MPI_INT, recv_rank, 0, MPI_COMM_WORLD); 
              MPI_Send(&num_to_next,  1, MPI_INT, recv_rank, 0, MPI_COMM_WORLD); 
                          
              int src_position = vertices_local_pL_offset[lvl+1] - 7 - num_to_next; // sum of vertices incl thsi level          
              // MPI_Ssend(&mask[src_position],  num_to_next, MPI_UINT8_T, recv_rank, 0, MPI_COMM_WORLD); 
              MPI_Send(&mask[src_position],  num_to_next, MPI_UINT8_T, recv_rank, 1, MPI_COMM_WORLD); 
            }     
            if (pid() == recv_rank) {
              MPI_Recv(&num_from_prev, 1, MPI_INT, send_rank, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); 
              MPI_Recv(&tmp[0], num_from_prev, MPI_UINT8_T, send_rank, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);           
            }
          }                        
        }
        
            
        if (vertices_local_corr %8 != 0) // sanity check
          MPI_Abort(MPI_COMM_WORLD, 2);
          
        //* TEMP EXPORT FILE STRUCTURE */
        #if 0
        FILE * fpt;
        char exp_name[80];
        sprintf(exp_name, "structure_pid%04i.dat", pid());
        fpt = fopen(exp_name, "w");
        fprintf(fpt, "%li\n",  length_w_spacing);
        for (int lvl = 0; lvl < grid->maxdepth; ++lvl) {       
          fprintf(fpt, "%i\n",  vertices_local_pL_offset_corr[lvl]);
          fprintf(fpt, "%i\n",  vertices_local_pL_corr[lvl]) ;
        }
        fclose(fpt);
        #endif
        //* END TEMP EXPORT FILE STRUCTURE */
        
        // zero copy bitmask creation
        
        int i = 0, cnt = 0;
        for (int lvl = 0; lvl < grid->maxdepth; ++lvl) {      
          int displacement = vertices_local_pL_offset_corr[lvl]; // offset of correct array within my array
          int count = vertices_local_pL_corr[lvl]; // corrected count        
          for (int c = 0; c < count; ++c) {
            mask[i] |= mask[displacement + c] << (7-cnt);        
            if (++cnt%8 == 0) {
               mask[++i] = 0; // incr i and clear next byte
               cnt=0; 
             }
          }
        }

    We need to calculate the global offset of the bit packets for the * MPI_Type_indexed datatype describing the tree data

        int vertices_global_offset_corr[grid->maxdepth]; // not "+1"
        vertices_global_offset_corr[0] = 0;
        unsigned int carryover = 0;
        for (int lvl = 0; lvl < grid->maxdepth; ++lvl){
        MPI_Exscan (&vertices_local_pL_corr[lvl], \
                    &vertices_global_offset_corr[lvl], \
                    1, \
                    MPI_INT, \
                    MPI_SUM, \
                    MPI_COMM_WORLD);    
          // Send Data to Process 0 for the Exscan on next Level
          if (pid() == (npe() - 1)) {
            unsigned int next_offset;
            next_offset = vertices_global_offset_corr[lvl] + vertices_local_pL_corr[lvl];
            MPI_Ssend(&next_offset, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);            
          }
          if (pid() == 0) {
            vertices_local_pL_corr[lvl] -= carryover; // remove temporary offset added in previous level
            
            MPI_Recv(&carryover, 1, MPI_INT, npe() - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);            
            if (lvl + 1 < grid->maxdepth ) {
              vertices_local_pL_corr[lvl+1] += carryover; // use array temporarily to carry over offset to next level
              vertices_global_offset_corr[lvl+1] = carryover;
            }      
          } 
        }

    Create a struct which will be written to RAM

        struct descBit_t{
          u_int32_t size;
          u_int8_t* data;
        }descBit_struct;
            
        descBit_struct.size = (descBits/8 + 1) * cell_size;
        descBit_struct.data =  &mask[0];

    Define where the data is located in Memory

        MPI_Aint m_displacements[2];
        MPI_Aint base_address;
        MPI_Get_address(&descBit_struct,         &base_address);
        MPI_Get_address(&descBit_struct.size,    &m_displacements[0]);
        MPI_Get_address(&descBit_struct.data[0], &m_displacements[1]);
        m_displacements[0] = MPI_Aint_diff(m_displacements[0], base_address);
        m_displacements[1] = MPI_Aint_diff(m_displacements[1], base_address);
        
        MPI_Datatype m_types[2] = { MPI_UINT32_T, MPI_BYTE };    
        //~

    Tree has to have some data, else there is a segmentationfault when setting view

        int m_lengths[2] = {1, vertices_local_corr/8}; // Bits
            
        MPI_Datatype m_view; // memory view
        MPI_Type_create_struct(2, m_lengths, m_displacements, m_types, &m_view);
        MPI_Type_commit(&m_view);    
        
        int lengths[grid->maxdepth];
        int displacements[grid->maxdepth];

    Define the File view for writing the bitmask [1/2]

        for(int lvl = 0; lvl < grid->maxdepth; ++lvl) {
          lengths[lvl]       = (int)vertices_local_pL_corr[lvl]/8; // bits, not bytes
          displacements[lvl] = (int)vertices_global_offset_corr[lvl]/8; // bits, not bytes
        }    
        
        MPI_Datatype tree_type_descBit;
        MPI_Type_indexed(grid->maxdepth, lengths, displacements, MPI_BYTE, &tree_type_descBit);
        MPI_Type_commit(&tree_type_descBit);

    Define the File view for writing the bitmask [2/2]

        MPI_Aint f_displacements[2] = {0, sizeof(u_int32_t)};
        int f_lengths[2] = {1, 1};
        MPI_Datatype f_types[2] = { MPI_UINT32_T, tree_type_descBit};    
        
        MPI_Datatype f_view; // file_view
        MPI_Type_create_struct(2, f_lengths, f_displacements, f_types, &f_view);
        MPI_Type_commit(&f_view);

    Set the view according to f_fiew (file_view)

        MPI_File_set_view(fp, offset, f_view, f_view, "native", MPI_INFO_NULL);

    Write from memory according to m_view (memory_view)

        MPI_File_write_all(fp, &descBit_struct, 1, m_view, MPI_STATUS_IGNORE);     
        
        offset += (descBits/8+1) * sizeof(u_int8_t) + sizeof(u_int32_t);
        
        MPI_Type_free(&m_view);
        MPI_Type_free(&f_view);           
        MPI_Type_free(&tree_type_descBit);
        
        free(mask); mask = NULL;       
        descBit_struct.data = NULL;   
      }

    Write LEVEL

    #if WRITE_HTG_LEVEL
      {
        struct level_t{
          u_int32_t size;
          u_int8_t* data;
        }level_struct;
            
        cell_size=sizeof(u_int8_t);
        level_struct.size =  vertices * cell_size;
        level_struct.data = (u_int8_t*) malloc(vertices_local*cell_size);
     
        MPI_Aint m_displacements[2];
        MPI_Aint base_address;
        MPI_Get_address(&level_struct,         &base_address);
        MPI_Get_address(&level_struct.size,    &m_displacements[0]);
        MPI_Get_address(&level_struct.data[0], &m_displacements[1]);
        m_displacements[0] = MPI_Aint_diff(m_displacements[0], base_address);
        m_displacements[1] = MPI_Aint_diff(m_displacements[1], base_address);
        
        MPI_Datatype m_types[2] = { MPI_UINT32_T, MPI_UINT8_T};
        int m_lengths[2] = {1, vertices_local};
        
        MPI_Datatype m_view;
        MPI_Type_create_struct(2, m_lengths, m_displacements, m_types, &m_view);
        MPI_Type_commit(&m_view);    
        
        int lengths[grid->maxdepth+1];
        int displacements[grid->maxdepth+1];
        for(int lvl = 0; lvl < grid->maxdepth + 1; ++lvl){
          lengths[lvl]       = (int) vertices_local_pL[lvl];
          displacements[lvl] = (int) vertices_global_offset[lvl];
        }
        MPI_Datatype tree_type_level;
        MPI_Type_indexed(grid->maxdepth + 1, lengths, displacements, MPI_UINT8_T, &tree_type_level);
        MPI_Type_commit(&tree_type_level);
        
        MPI_Aint f_displacements[2] = {0, sizeof(u_int32_t)};
        int f_lengths[2] = {1, 1};
        MPI_Datatype f_types[2] = { MPI_UINT32_T, tree_type_level};
        
        MPI_Datatype f_view;
        MPI_Type_create_struct(2, f_lengths, f_displacements, f_types, &f_view);
        MPI_Type_commit(&f_view);
        
        long index = 0;   
        for(int lvl=0; lvl<=grid->maxdepth;++lvl)		
          foreach_level(lvl,serial,noauto) 
            if (is_local(cell)) 
              level_struct.data[index++] = (u_int8_t)lvl;
                              
        MPI_File_set_view(fp, offset, f_view,f_view, "native", MPI_INFO_NULL);        
        MPI_File_write_all(fp, &level_struct, 1, m_view, MPI_STATUS_IGNORE);
    
        offset += vertices*cell_size + sizeof(u_int32_t);
    
        free(level_struct.data); level_struct.data = NULL;
        MPI_Type_free(&m_view);       
        MPI_Type_free(&f_view);       
        MPI_Type_free(&tree_type_level);
      }
    #endif // end WRITE_HTG_LEVEL
      {
        struct scalar_t{
          u_int32_t size;
          float* data;
        }scalar_struct;
    
        cell_size=sizeof(float);
        scalar_struct.size = vertices*cell_size;
        scalar_struct.data = (float*) malloc(vertices_local*cell_size);
        
        MPI_Aint m_displacements[2];
        MPI_Aint base_address;
        MPI_Get_address(&scalar_struct,         &base_address);
        MPI_Get_address(&scalar_struct.size,    &m_displacements[0]);
        MPI_Get_address(&scalar_struct.data[0], &m_displacements[1]);
        m_displacements[0] = MPI_Aint_diff(m_displacements[0], base_address);
        m_displacements[1] = MPI_Aint_diff(m_displacements[1], base_address);
        
        MPI_Datatype m_types[2] = { MPI_UINT32_T, MPI_FLOAT};
        int m_lengths[2] = {1, vertices_local};
        
        MPI_Datatype m_view;
        MPI_Type_create_struct(2, m_lengths, m_displacements, m_types, &m_view);
        MPI_Type_commit(&m_view);    
        
        int lengths[grid->maxdepth+1];
        int displacements[grid->maxdepth+1];
        for(int lvl = 0; lvl < grid->maxdepth + 1; ++lvl){
          lengths[lvl]       = (int) vertices_local_pL[lvl];
          displacements[lvl] = (int) vertices_global_offset[lvl];
        }
        MPI_Datatype tree_type_scalar;
        MPI_Type_indexed(grid->maxdepth + 1, lengths, displacements, MPI_FLOAT, &tree_type_scalar);
        MPI_Type_commit(&tree_type_scalar);
        
        MPI_Aint f_displacements[2] = {0, sizeof(u_int32_t)};
        int f_lengths[2] = {1, 1};
        MPI_Datatype f_types[2] = { MPI_UINT32_T, tree_type_scalar};
        
        MPI_Datatype f_view;
        MPI_Type_create_struct(2, f_lengths, f_displacements, f_types, &f_view);
        MPI_Type_commit(&f_view);
        
        for (scalar s in list) {
          long index = 0;
          for(int lvl=0; lvl<=grid->maxdepth;++lvl)
            foreach_level(lvl,serial,noauto)
              if (is_local(cell))
                 scalar_struct.data[index++] = (float)val(s);
          
          MPI_File_set_view(fp, offset, f_view,f_view, "native", MPI_INFO_NULL);      
          MPI_File_write_all(fp, &scalar_struct, 1, m_view, MPI_STATUS_IGNORE); 
          offset += vertices*cell_size + sizeof(u_int32_t);
        }
        free(scalar_struct.data); scalar_struct.data = NULL;
        MPI_Type_free(&m_view);       
        MPI_Type_free(&f_view);
        MPI_Type_free(&tree_type_scalar);  
      }   
      {
        struct vector_t{
          u_int32_t size;
          float* data;
        }vector_struct;
    
        cell_size=3*sizeof(float);
        vector_struct.size = vertices*cell_size;
        vector_struct.data = (float*) malloc(vertices_local*cell_size);
        
        MPI_Aint m_displacements[2];
        MPI_Aint base_address;
        MPI_Get_address(&vector_struct,         &base_address);
        MPI_Get_address(&vector_struct.size,    &m_displacements[0]);
        MPI_Get_address(&vector_struct.data[0], &m_displacements[1]);
        m_displacements[0] = MPI_Aint_diff(m_displacements[0], base_address);
        m_displacements[1] = MPI_Aint_diff(m_displacements[1], base_address);
        
        MPI_Datatype m_types[2] = { MPI_UINT32_T, MPI_FLOAT};
        int m_lengths[2] = {1, 3*vertices_local};
        
        MPI_Datatype m_view;
        MPI_Type_create_struct(2, m_lengths, m_displacements, m_types, &m_view);
        MPI_Type_commit(&m_view);    
        
        int lengths[grid->maxdepth+1];
        int displacements[grid->maxdepth+1];
        for(int lvl = 0; lvl < grid->maxdepth + 1; ++lvl) {
          lengths[lvl]       = 3*(int)vertices_local_pL[lvl];
          displacements[lvl] = 3*(int)vertices_global_offset[lvl];
        }
        MPI_Datatype tree_type_vector;
        MPI_Type_indexed(grid->maxdepth + 1, lengths, displacements, MPI_FLOAT, &tree_type_vector);
        MPI_Type_commit(&tree_type_vector);
        
        MPI_Aint f_displacements[2] = {0, sizeof(u_int32_t)};
        int f_lengths[2] = {1, 1};
        MPI_Datatype f_types[2] = { MPI_UINT32_T, tree_type_vector};
        
        MPI_Datatype f_view;
        MPI_Type_create_struct(2, f_lengths, f_displacements, f_types, &f_view);
        MPI_Type_commit(&f_view);        
    
        for (vector v in vlist) {
          long index = 0;
            for(int lvl=0; lvl<=grid->maxdepth;++lvl)
              foreach_level(lvl,serial,noauto) 
                if (is_local(cell)) {
                  #if dimension == 2					
                  vector_struct.data[index]   = (float)val(v.y);
                  vector_struct.data[index+1] = (float)val(v.x);
                  vector_struct.data[index+2] = (float)0.;
                  index += 3;
                  #elif dimension == 3					
                  vector_struct.data[index]   = (float)val(v.z);
                  vector_struct.data[index+1] = (float)val(v.y);
                  vector_struct.data[index+2] = (float)val(v.x);
                  index += 3;
                  #endif					
                } 			            
          MPI_File_set_view(fp, offset, f_view,f_view, "native", MPI_INFO_NULL);
          MPI_File_write_all(fp, &vector_struct, 1, m_view, MPI_STATUS_IGNORE);       
          offset += vertices*cell_size + sizeof(u_int32_t);
        }
        free(vector_struct.data); vector_struct.data = NULL;
        MPI_Type_free(&m_view);       
        MPI_Type_free(&f_view);
        MPI_Type_free(&tree_type_vector);
      } 
      
      MPI_File_set_view(fp, offset, MPI_BYTE,MPI_BYTE, "native", MPI_INFO_NULL);
    	
    	if (pid() == 0)	
    		Write2File(sprintf(buffer, "ENDBINARY\n\t</AppendedData>\n</VTKFile>\n"));    
      
    #if HTG_SPEED_STATS  
      MPI_Barrier(MPI_COMM_WORLD);  
    	double end = MPI_Wtime();  
    	fprintf (stdout,"Write Time Taken: %f ms, Throughput: %.2f MB/s, Stats %.2f%%, Header: %.2f%%\n",\
        (end-start)*1000., (double)(offset+strlen(buffer))/(end-start)/(double)sq(1024),\
        t_stats/(end-start)*100.,t_header/(end-start)*100.);    
      fflush(stdout);
    #endif    
      MPI_File_sync(fp);
      MPI_Barrier(MPI_COMM_WORLD);  
      #if defined(_OPENMP)
        omp_set_num_threads(num_omp);
      #endif
    }
    #else
    void output_htg_data(scalar * list, vector * vlist, FILE * fp)
    {
    	#if defined(_OPENMP)
    		int num_omp = omp_get_max_threads();
    		omp_set_num_threads(1);
    	#endif  
    	
      unsigned int vertices_local = 0;
    	unsigned int descBits_local = 0;  
      unsigned int vertices_local_pL[grid->maxdepth+1];  // pL = per Level 
    
        
    	for (int lvl = 0; lvl < grid->maxdepth+1; ++lvl){
        vertices_local_pL[lvl] = 0;
    		foreach_level(lvl,serial,noauto) 
    			if(is_local(cell)) 
    				vertices_local_pL[lvl]++;
        
        vertices_local += vertices_local_pL[lvl];         
      }
      
      descBits_local = vertices_local - vertices_local_pL[grid->maxdepth];
      
    #if HTG_SPEED_STATS
      struct timespec start;
      
      if( clock_gettime( CLOCK_REALTIME, &start) == -1 ) {
        perror( "clock gettime" );    
      }
    #endif
    • Calculate min,max for xml header (prob. not needed?)
    #if HTG_SPEED_STATS       
      struct timespec t_stats_start;  
      if( clock_gettime( CLOCK_REALTIME, &t_stats_start) == -1 ) {
        perror( "clock gettime" );
         
      } 
    #endif  
      double min_val[list_len(list)];
      double max_val[list_len(list)];
      {
        int i = 0;
        for (scalar s in list) {
    #if HEADER_MIN_MAX_VAL
          stats stat = statsf(s);
          min_val[i]=stat.min;
          max_val[i]=stat.max;
    #else
          min_val[i]=0.;
          max_val[i]=0.;
    #endif
          i++;
        }
    	}
      double min_val_v[vectors_len(vlist)];
      double max_val_v[vectors_len(vlist)];
      { 
        int i = 0;	
        for (vector v in vlist) {
    #if HEADER_MIN_MAX_VAL
          min_val_v[i]= HUGE;
          max_val_v[i]= -HUGE;
          foreach_dimension(){
            stats stat = statsf(v.x);	
            min_val_v[i] = min(stat.min,min_val_v[i]);
            max_val_v[i]=  max(stat.max,max_val_v[i]);			
          }
    #else
          min_val_v[i]= 0.;
          max_val_v[i]= 0.;
    #endif
          i++;
        }
    	}
    #if HTG_SPEED_STATS    
      struct timespec t_stats_stop;  
      if( clock_gettime( CLOCK_REALTIME, &t_stats_stop) == -1 ) {
        perror( "clock gettime" );
         
      }  
      // to do t_stats_stop    
      double t_stats = ( t_stats_stop.tv_sec - t_stats_start.tv_sec )
                      + (double)( t_stats_stop.tv_nsec - t_stats_start.tv_nsec )
                        / (double)1000000000L;
    	
      struct timespec t_header_start;
      if( clock_gettime( CLOCK_REALTIME, &t_header_start) == -1 ) {
        perror( "clock gettime" );
         
      }  
    #endif

    File Header vtkHyperTreeGrid

        //~ WriteHeader(&offset, infoStruct);
      int maj_v = VTK_FILE_VERSION/10, min_v = VTK_FILE_VERSION%10;
      
      fprintf(fp,\
        "<VTKFile %s version=\"%i.%i\" %s %s>\n",\
          "type=\"HyperTreeGrid\"", \
          maj_v, min_v , \
          "byte_order=\"LittleEndian\" ",\
          "header_type=\"UInt32\""); 
    
    #if dimension==2
      fprintf(fp,\
        "\t<HyperTreeGrid BranchFactor=\"2\" TransposedRootIndexing=\"0\" Dimensions=\"%d %d %d\">\n",\
        2,2,1);    
    #elif dimension==3
      fprintf(fp,\
        "\t<HyperTreeGrid BranchFactor=\"2\" TransposedRootIndexing=\"0\" Dimensions=\"%d %d %d\">\n", 2,2,2);
    #endif
      fprintf(fp,"\t\t<Grid>\n");
    #if dimension==2
      fprintf(fp, "\t\t\t<DataArray type=\"Float64\" Name=\"XCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",Y0, Y0+L0);
      fprintf(fp, "\t\t\t\t%g %g",Y0, Y0+L0);
      fprintf(fp,"\n\t\t\t</DataArray>\n");
      fprintf(fp, "\t\t\t<DataArray type=\"Float64\" Name=\"YCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",X0, X0+L0);
      fprintf(fp, "\t\t\t\t%g %g",X0, X0+L0);
      fprintf(fp,"\n\t\t\t</DataArray>\n");	  
      fprintf(fp, "\t\t\t<DataArray type=\"Float64\" Name=\"ZCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",Z0, Z0+L0);		
      fprintf(fp, "\t\t\t\t%g %g", 0., 0.);	  
    #elif dimension==3
      fprintf(fp, "\t\t\t<DataArray type=\"Float64\" Name=\"XCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",\
          Z0, Z0+L0);
      fprintf(fp, "\t\t\t\t%g %g",Z0, Z0+L0);
      fprintf(fp,"\n\t\t\t</DataArray>\n");
      fprintf(fp, "\t\t\t<DataArray type=\"Float64\" Name=\"YCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",\
          Y0, Y0+L0);
      fprintf(fp, "\t\t\t\t%g %g",Y0, Y0+L0);
      fprintf(fp,"\n\t\t\t</DataArray>\n");	  
      fprintf(fp, "\t\t\t<DataArray type=\"Float64\" Name=\"ZCoordinates\" NumberOfTuples=\"2\" format=\"ascii\" RangeMin=\"%g\" RangeMax=\"%g\">\n",\
          X0, X0+L0);		
      fprintf(fp, "\t\t\t\t%g %g",X0, X0+L0);	  
    #endif
      fprintf(fp,"\n\t\t\t</DataArray>\n");	  
      fprintf(fp,"\t\t</Grid>\n");
      fprintf(fp,"\t\t<Trees>\n");
      
      // Tree Begin
      unsigned int byte_offset = 0;
    #if VTK_FILE_VERSION == 10
      fprintf(fp,"\t\t\t<Tree Index=\"0\" NumberOfLevels=\"%d\" NumberOfVertices=\"%u\">\n",\
          grid->maxdepth + 1, vertices_local);   	
      
      // Array Descriptor Bits    
      fprintf(fp,"\t\t\t\t<DataArray type=\"Bit\" Name=\"Descriptor\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"0\" RangeMax=\"1\" offset=\"%u\"/>\n",\
        descBits_local, byte_offset);        
      byte_offset += (descBits_local/8+1) * sizeof(u_int8_t) + sizeof(u_int32_t);
      
      // Array NbVerticesByLevel
      fprintf(fp,"\t\t\t\t<DataArray type=\"Int64\" Name=\"NbVerticesByLevel\" NumberOfTuples=\"%d\" format=\"ascii\" RangeMin=\"1\" RangeMax=\"%u\" >\n\t\t\t\t\t",\
          grid->maxdepth + 1, vertices_local_pL[grid->maxdepth]);
          
      for (int lvl = 0; lvl <= grid->maxdepth; lvl++) {
        fprintf(fp,"%u ", vertices_local_pL[lvl]);
      }		
      fprintf(fp,"\n\t\t\t\t</DataArray>\n");
    #endif     
    
    #if VTK_FILE_VERSION == 20
      // Array Descriptor Bits    
      fprintf(fp,"\t\t\t\t<DataArray type=\"Bit\" Name=\"Descriptors\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"0\" RangeMax=\"1\" offset=\"%u\"/>\n",\
        descBits_local, byte_offset);        
      byte_offset += (descBits_local/8+1) * sizeof(u_int8_t) + sizeof(u_int32_t);
      
      fprintf(fp,"\t\t\t\t<DataArray type=\"Int64\" Name=\"NumberOfVerticesPerDepth\" NumberOfTuples=\"%d\" format=\"ascii\" RangeMin=\"1\" RangeMax=\"%u\" >\n\t\t\t\t\t",\
          grid->maxdepth + 1, vertices_local_pL[grid->maxdepth]);
            
      for (int lvl = 0; lvl <= grid->maxdepth; lvl++) {
        fprintf(fp,"%u ", vertices_local_pL[lvl]);
      }		
      fprintf(fp,"\n\t\t\t\t</DataArray>\n");
      
      fprintf(fp,"\t\t\t\t<DataArray type=\"Int64\" Name=\"TreeIds\" NumberOfTuples=\"1\" format=\"ascii\" RangeMin=\"0\" RangeMax=\"0\" >\n");
      fprintf(fp,"\t\t\t\t\t%i\n",0);
      fprintf(fp,"\t\t\t\t</DataArray>\n");
      
      fprintf(fp,"\t\t\t\t<DataArray type=\"UInt32\" Name=\"DepthPerTree\" NumberOfTuples=\"1\" format=\"ascii\" RangeMin=\"%i\" RangeMax=\"%i\" >\n",\
      grid->maxdepth+1, grid->maxdepth+1);
      fprintf(fp,"\t\t\t\t\t%i\n",grid->maxdepth+1);
      fprintf(fp,"\t\t\t\t</DataArray>\n");
              
      fprintf(fp,"\t\t</Trees>\n");
    #endif  
      // Cell Data Begin
      fprintf(fp,"\t\t\t\t<CellData>\n");
      
    #if WRITE_HTG_LEVEL    
      fprintf(fp,"\t\t\t\t\t<DataArray type=\"UInt8\" Name=\"Level\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"0\" RangeMax=\"%i\" offset=\"%u\"/>\n",\
          vertices_local, grid->maxdepth, byte_offset);		
      byte_offset += vertices_local * sizeof(u_int8_t) + sizeof(u_int32_t);
    #endif    
      {
        int i = 0;
        for (scalar s in list)
        {			
          fprintf(fp,"\t\t\t\t\t<DataArray type=\"Float32\" Name=\"%s\" NumberOfTuples=\"%u\" format=\"appended\" RangeMin=\"%g\" RangeMax=\"%g\" offset=\"%u\"/>\n",\
              s.name, vertices_local, min_val[i], max_val[i], byte_offset);			    
          byte_offset += vertices_local * sizeof(float_t) + sizeof(u_int32_t);
          i++;
        }
      }
      {
        int i = 0;
        for (vector v in vlist)
        {
          char *vname = strtok(v.x.name, ".");
          fprintf(fp, "\t\t\t\t\t<DataArray type=\"Float32\" NumberOfComponents=\"%i\" Name=\"%s\" NumberOfTuples=\"%u\" format=\"appended\"  RangeMin=\"%g\" RangeMax=\"%g\" offset=\"%u\"/>\n",\
              3,vname, vertices_local, min_val_v[i],max_val_v[i],byte_offset);
          byte_offset += vertices_local * 3 * sizeof(float_t) + sizeof(u_int32_t);
          i++;
        } 
      }  
      fprintf(fp,"\t\t\t\t</CellData>\n");
    #if VTK_FILE_VERSION == 10
      fprintf(fp,"\t\t\t</Tree>\n\t\t</Trees>\n");
    #endif
    
      fprintf(fp,"\t</HyperTreeGrid>\n\t<AppendedData encoding=\"raw\">\n_");
    
    #if HTG_SPEED_STATS    
      struct timespec t_header_stop;
      if( clock_gettime( CLOCK_REALTIME, &t_header_stop) == -1 ) {
        perror( "clock gettime" );
      } 
      
      double t_header = ( t_header_stop.tv_sec - t_header_start.tv_sec )
                      + (double)( t_header_stop.tv_nsec - t_header_start.tv_nsec )
                        / (double)1000000000L;
    #endif

    Write Bitmask

      int cell_size;
    
    	cell_size=sizeof(u_int8_t);    
    	  
      int vertices_local_corr = ((descBits_local/8)+1)*8;
    
      u_int32_t prepend_size = vertices_local_corr;
      fwrite(&prepend_size, sizeof(u_int32_t), 1, fp); 		
      u_int8_t* write_cache = (u_int8_t*)calloc(vertices_local_corr,cell_size);
      long index = 1;
      for(int lvl=0; lvl<grid->maxdepth;++lvl){ // Bitmask ist "0" auf grid->maxdepth	und muss nicht geschrieben werden
    		foreach_level(lvl,serial,noauto)
    			if (is_local(cell)) {
    				if(is_leaf(cell)){					
    					write_cache[index++] = 0; // ascii: 0					
    				} else {					
    					write_cache[index++] = 1; // ascii: 1
    				}
          }
    	}
    
      // Byte to bits
      for (int i = 0; i < vertices_local_corr/8 ; ++i){
        for ( int j = 0; j < 8; ++j ) {
          write_cache[i] |= write_cache[(8*i+j) + 1] << (7-j);        
          if ((j+1) == 8) {
             write_cache[i+1] = 0; // clear next byte         
          }
        }
      }
      
      fwrite(&write_cache[0], cell_size, vertices_local_corr/8, fp);
      free(write_cache);
      write_cache = NULL;

    Write LEVEL

    #if WRITE_HTG_LEVEL
     
     	cell_size=sizeof(u_int8_t);
    	
      prepend_size = vertices_local * cell_size;
      fwrite(&prepend_size, sizeof(u_int32_t), 1, fp); 		
    
    	for(int lvl=0; lvl<=grid->maxdepth;++lvl){
    		u_int8_t * write_cache = (u_int8_t*) malloc(vertices_local_pL[lvl]*sizeof(u_int8_t));
    		long index = 0;
    		foreach_level(lvl,serial,noauto)
    			if (is_local(cell))
    				write_cache[index++] = (u_int8_t)lvl;
    
    		fwrite(&write_cache[0], cell_size, vertices_local_pL[lvl], fp);
    		free(write_cache);
    		write_cache = NULL;
    	}
    	
    #endif // end WRITE_HTG_LEVEL
      
    	for (scalar s in list) {
    		cell_size=sizeof(float_t);
    		  
        u_int32_t prepend_size = vertices_local * cell_size;
        fwrite(&prepend_size, sizeof(u_int32_t), 1, fp); 		  
    	  	
    		for(int lvl=0; lvl < grid->maxdepth + 1; ++lvl){
          
    			// offset per level
    			float_t * write_cache = (float_t*) malloc(vertices_local_pL[lvl]*cell_size);
    			long index = 0;
          
    			foreach_level(lvl,serial,noauto)
    				if (is_local(cell)) 			
    					write_cache[index++] = val(s);
    			
    			fwrite(&write_cache[0], cell_size, vertices_local_pL[lvl], fp);      
    			free(write_cache);
    			write_cache = NULL;	        
        }
      }
      
    	for (vector v in vlist) {
    		cell_size = 3*sizeof(float_t);
    		
        u_int32_t prepend_size = vertices_local * cell_size;
        fwrite(&prepend_size, sizeof(u_int32_t), 1, fp); 		  
    			
    		for(int lvl=0; lvl<=grid->maxdepth;++lvl){
    			
    			float_t * write_cache = (float_t*) malloc(vertices_local_pL[lvl]*cell_size);
    			long index = 0;
    			foreach_level(lvl,serial,noauto)
    				if (is_local(cell)) {
    					#if dimension == 2
    			
    					write_cache[index] = val(v.y);
    					write_cache[index+1] = val(v.x);
    					write_cache[index+2] = 0.;
    					index += 3;
    					#elif dimension == 3
    			
    					write_cache[index] = val(v.z);
    					write_cache[index+1] = val(v.y);
    					write_cache[index+2] = val(v.x);
    					index += 3;
    					#endif	
            }
    			fwrite(&write_cache[0], cell_size, vertices_local_pL[lvl], fp);
    			free(write_cache);
    			write_cache = NULL;	
    		}    	
    	} 
    
    	fprintf(fp, "ENDBINARY\n\t</AppendedData>\n</VTKFile>\n");    
      
    #if HTG_SPEED_STATS  
     struct timespec stop;
     if( clock_gettime( CLOCK_REALTIME, &stop) == -1 ) {
        perror( "clock gettime" );    
      }
      
      long offset = ftell(fp);
      
      double t_file = ( stop.tv_sec - start.tv_sec )
                      + (double)( stop.tv_nsec - start.tv_nsec )
                        / (double)1000000000L;
    	fprintf (stdout,"Write Time Taken: %lf ms, Throughput: %.2f MB/s, Stats %.2f%%, Header: %.2f%%\n",\
        t_file*1000., (double)offset/t_file/(double)sq(1024),\
        t_stats/t_file*100.,t_header/t_file*100.);
      fflush(stdout);
    #endif
      fflush(fp);
      #if defined(_OPENMP)
        omp_set_num_threads(num_omp);
      #endif
      
    }
    #endif