r2mfile.c

用于2维的射线追踪

				   ci->item[i].in->prof,
				   ci->item[i].out->lat,
				   ci->item[i].out->lon,
				   ci->item[i].out->prof,
				   ci->item[i].P_length,
				   ci->item[i].S_length);
		    }
		    my_fprintf(fd, "\n");
		    mesh_cellinfo_free(c[z][x][y]);
		    c[z][x][y] = NULL;
		}
	    }
	}
    }
    my_fprintf(fd, "# end\n");
    my_fclose(fd);
#ifdef DEBUG
    date_stamp = get_date_stamp();
    fprintf(stdout, "[%s] make_domain_info_file() ends\n", date_stamp);
    fflush(stdout);
    free(date_stamp);
#endif

    return (nb_cell_cpt);
}



/** \brief read a integer into the buffer **/
void read_int(int *val)
{
    char *buf;

    buf = strtok(NULL, " \t\n");
    *val = atoi(buf);
}

/** \brief read a double into the buffer **/
void read_double(double *val)
{
    char *buf;

    buf = strtok(NULL, " \t\n");
    *val = atof(buf);
}

/** \brief read a long into the buffer **/
void read_long(long int *val)
{
    char *buf;

    buf = strtok(NULL, " \t\n");
    *val = atol(buf);
}

/** \brief read a string into the buffer **/
char *read_string(void)
{
    return (strtok(NULL, " \t\n"));
}

/** \brief go to next line **/
void pass_line(char *buf)
{
    strtok(buf, "\n");
}
/** \brief import a fomated buffer into one cell.
 *
 * Reads cell data from a r2m text. Processes the
 * contents of one cell : fields face hit, block hits, items with rays in,out, 
 * P and S lengths, and ray id.
 *
 * @param c an allocated cell to store data into  
 * @param buff the buffer with the r2m data
 */
void import_cell_data_buffer( struct cell_info_t *c, char *buff)		
{
struct coord_geo_t in, out;
double Pl, Sl;
int nb_block_hit, nb_item;
int face_hit, block_hit;
int i;
long int rayid;


	/* face hit */
	for (i = 0; i < 6; i++) {
	    read_int(&face_hit);
	    c->faces_hit[i] += face_hit;
	}
	/* block hit */
	read_int(&nb_block_hit);
	for (i = 0; i < nb_block_hit; i++) {
	    read_int(&block_hit);
	    mesh_cellinfo_block_feed(c, block_hit);
	}

	/* item */
	read_int(&nb_item);
	for (i = 0; i < nb_item; i++) {
	    read_long(&rayid);	/* ray Id  */
	    read_double(&(in.lat));	/* lat in  */
	    read_double(&(in.lon));	/* lon in  */
	    read_double(&(in.prof));	/* prf in  */
	    read_double(&(out.lat));	/* lat out */
	    read_double(&(out.lon));	/* lon out */
	    read_double(&(out.prof));	/* prf out */
	    read_double(&Pl);	/* P length  */
	    read_double(&Sl);	/* S length  */

	    /* add the item in the current cell */
	    mesh_cellinfo_length_feed(c, &in, &out, Pl, Sl, rayid);
      }
}

/** \brief Only read data into the buffer, but don't do anything else */
void fake_import_cell_data_buffer( struct cell_info_t *c, char *buff)		
{
struct coord_geo_t in, out;
double Pl, Sl;
int nb_block_hit, nb_item;
int face_hit, block_hit;
int i;
long int rayid;


	/* face hit */
	for (i = 0; i < 6; i++) {
	    read_int(&face_hit);
	}
	/* block hit */
	read_int(&nb_block_hit);
	for (i = 0; i < nb_block_hit; i++) {
	    read_int(&block_hit);
	}

	/* item */
	read_int(&nb_item);
	for (i = 0; i < nb_item; i++) {
	    read_long(&rayid);	/* ray Id  */
	    read_double(&(in.lat));	/* lat in  */
	    read_double(&(in.lon));	/* lon in  */
	    read_double(&(in.prof));	/* prf in  */
	    read_double(&(out.lat));	/* lat out */
	    read_double(&(out.lon));	/* lon out */
	    read_double(&(out.prof));	/* prf out */
	    read_double(&Pl);	/* P length  */
	    read_double(&Sl);	/* S length  */
      }
}


/** \brief import a r2m fomated buffer into a light mesh.
 *
 * Used to read a r2m formated text, i.e. with the number of cells 
 * preceding the cells data blocks. 
 * For each cell data block, allocates a cell in a light mesh,
 * and call the auxiliary function import_cell_data_buffer()
 * to read and store the data.
 *
 * Returns the number of cell read.
 *
 @param c the light mesh
 @param buff the buffer with the r2m data
 @param size buffer size
 @param first_layer filter layers
 @param last_layer filter layers
 @param rank the current process id (for trace purpose only)
 @param dom_id the current domain id (for trace purpose only)
 */
int import_cell_buffer(struct cell_info_t ****c,	
			     char *buff,		
			     int size,	
			     int first_layer,
			     int last_layer,
			     int rank,
			     int dom_id)
{
    int nb_cell;
    int n, x, y, z;
#ifdef DEBUG
    char *date_stamp;
#endif 

    if ((buff == NULL) || (size==0)) {
#ifdef DEBUG
    fprintf(stderr,
         "[process %d] %s:import_cell_buffer():%d dom_id=%d : empty buffer !\n",
	    rank,__FILE__,__LINE__,dom_id);
#endif 
	return(0);
    }
    pass_line(buff);
    read_int(&nb_cell);
    /*assert(nb_cell>0);*/
    
#ifdef DEBUG
    date_stamp = get_date_stamp();
    fprintf(stderr,
	    "[process %d] %s:import_cell_buffer():%d  dom_id=%d, nb_cell=%d *** (%s)\n",
	    rank, __FILE__,__LINE__,dom_id, nb_cell,date_stamp);
    free(date_stamp);
#endif

    for (n = 0; n < nb_cell; n++) {
	/* cell coord */
	read_int(&z);
	read_int(&x);
	read_int(&y);

	/* filter layer */
	if ((first_layer >= 0) && (last_layer >= 0)) {
	    if ((z < first_layer) || (z > last_layer)) {
		fake_import_cell_data_buffer(c[z][x][y], buff);
		continue;
	    }
	}

	if (!c[z][x][y]) 
	    c[z][x][y] = mesh_cellinfo_onecell_alloc();

      import_cell_data_buffer( c[z][x][y] , buff);
    }

#ifdef DEBUG
    date_stamp = get_date_stamp();
    fprintf(stderr,
	    "[process %d] %s:import_cell_buffer():%d  dom_id=%d END *** (%s)\n",
	    rank, __FILE__,__LINE__,dom_id, date_stamp);
    free(date_stamp);
#endif
    return (nb_cell);
}


/** \brief merge two cells and returns the resulting cell.
 *  
 * The merge of two cells is an associative, commutative operation.
 * The function does not free the two cells passed as arguments.
 * @param c1 first cell to merge
 * @param c2 second cell to merge
 */
cell_info_t *merge_cells(cell_info_t *c1,cell_info_t *c2) {
cell_info_t *c;
int i;
	
	c = mesh_cellinfo_onecell_alloc();
      /* face hit */
	for (i = 0; i < 6; i++) {
	    c->faces_hit[i] = c1->faces_hit[i] + c2->faces_hit[i]; 
	}
	/* block hit */
	for (i = 0; i < c1->nblocks; i++) 
	    mesh_cellinfo_block_feed(c, c1->block_hit[i]);
	for (i = 0; i < c2->nblocks; i++) 
	    mesh_cellinfo_block_feed(c, c2->block_hit[i]);
      /* items */
      for (i = 0; i < c1->nitems; i++) {
	    mesh_cellinfo_length_feed (c, 
                                    c1->item[i].in, 
                                    c1->item[i].out, 
                                    c1->item[i].P_length,
                                    c1->item[i].S_length,
                                    c1->item[i].rayid);
      }
      for (i = 0; i < c2->nitems; i++) {
          mesh_cellinfo_length_feed(c, 
                                    c2->item[i].in, 
                                    c2->item[i].out, 
                                    c2->item[i].P_length,
                                    c2->item[i].S_length,
                                    c2->item[i].rayid);
      }
          
      return(c);
}

/*------------------------ZLIB----------------------------------------------*/
#ifdef USE_ZLIB
/** \brief import a compressed buffer and return a pointer to a uncompressed buffer  
 *
 * 
 * @param c the light mesh
 * @param zbuff the compressed buffer 
 * @param size compressed buffer size
 * @param uncompr_size uncompressed buffer size (returned)
 * @param rank the current process id 
 */
char *decompress_cell_data_gzbuffer(struct cell_info_t ****c,	
				    char *zbuff,		
				    int size,	
				    int *uncompr_size,	
				    int rank)
{
    int err;
    z_stream d_stream;		/* decompression stream */
    Byte *compr;
    Byte *uncompr;
    uLong comprLen;
    int count = 0;
#ifdef ZDEBUG
    FILE *fd;
    char filename[256];
#endif

    if (size == 0)
	return (NULL);

    /* Swap byte if necessary */
    memcpy((void *) uncompr_size, (void *) (zbuff + size - 4), 4);
    if (TestByteOrder() == MY_LITTLE_ENDIAN) {
	fprintf(stderr, "P[%d] Byte order is LITTLE_ENDIAN\n", rank);
    } else {
	fprintf(stderr, "P[%d] Byte order is BIG_ENDIAN\n", rank);
	*uncompr_size = SwapFourBytes(*uncompr_size);
    }


    fprintf(stdout,
	    "p[%d] data size uncompressed/compressed = %d/%d (%.1f)\n",
	    rank, *uncompr_size, size,
	    (float) (*uncompr_size) / (float) size);

    compr = (Byte *) (zbuff + 10);	/* skip the gzip header */
    size = -10;
    comprLen = (uLong) size;

    uncompr = (Byte *) calloc((uInt) (*uncompr_size + 1), 1);
    if (uncompr == NULL) {
#ifdef USE_MPI
	MPI_Abort(MPI_COMM_WORLD, 1);
#endif
	exit(1);
    }

    d_stream.zalloc = (alloc_func) 0;
    d_stream.zfree = (free_func) 0;
    d_stream.opaque = (voidpf) 0;

    d_stream.next_in = compr;
    d_stream.avail_in = 0;

    err = inflateInit2(&d_stream, -15);
    CHECK_ERR(err, "import_cell_buffer2: inflateInit");

#ifdef ZDEBUG
    snprintf(filename, 255, "gzbuffer-%.2d", rank);
    fd = fopen(filename, "w");
#endif

    while (d_stream.total_in < comprLen) {
	count++;
	d_stream.avail_in = d_stream.avail_out = *uncompr_size;
	d_stream.next_out = uncompr;
	err = inflate(&d_stream, /*Z_SYNC_FLUSH */ Z_FULL_FLUSH);
	uncompr[*uncompr_size] = '\0';
#ifdef ZDEBUG
	fprintf(fd, "%s", (char *) uncompr);
#endif
	if (err == Z_STREAM_END)
	    break;
	CHECK_ERR(err, "inflate");
    }

    err = inflateEnd(&d_stream);
    CHECK_ERR(err, "import_cell_buffer2: inflateEnd");
#ifdef ZDEBUG
    fclose(fd);
    fprintf(stdout, "p[%d] gzbuffer uncompressed\n", rank);
    fflush(stdout);
#endif
    return ((char *) uncompr);
}
#endif