r2mfile.c

用于2维的射线追踪

/*
 * Ray2mesh : software for geophysicists.
 * Compute various scores attached to the mesh cells, based on geometric
   information that rays bring when the traverse the cell.
 *
 * Copyright (C) 2003, St閜hane Genaud and Marc Grunberg
 *
 * This tool is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <sys/mman.h>

#ifdef USE_MPI
#include <mpi.h>
#endif

#ifdef USE_ZLIB
#include <zlib.h>
#endif

#include <mesh/mesh.h>
#include <mesh/cellinfo.h>
#include <mesh/layer.h>

#include "ray2mesh.h"
#include "util.h"
#include "r2mfile.h"
#include "merge.h"
#include "byte_order.h"

#ifdef USE_ZLIB
#define my_fprintf gzprintf
#define my_File    gzFile
#define my_fopen   gzopen
#define my_fclose  gzclose
#define openmode   "wb6f"
#else
#define my_fprintf fprintf
#define my_fopen   fopen
#define my_File    FILE
#define my_fclose  fclose
#define openmode   "w"
#endif

#ifdef USE_ZLIB
#define CHECK_ERR(err, msg) { \
    if (err != Z_OK) { \
        fprintf(stderr, "%s error: %d\n", msg, err); \
	MPI_Abort(MPI_COMM_WORLD, 1); \
        exit(1); \
    } \
}
#endif

/*#define ZDEBUG */

/** \brief return where starts and ends the domain "id" in layer "id" 
 *
 *@param m the mesh description 
 *@param layer_id the layer id
 *@param domain_id the domain id
 *@param nb_domain total number of domains 
 *@param from where the domain id starts
 *@param to where the domain id ends
 **/
void get_domain_bounds_bis (struct mesh_t *m,	
		       int layer_id,		
		       int domain_id,		
		       int nb_domain,		
		       int *from,		
		       int *to)
{
    /** where ends the domain id */
    int nb_lon_cell;
    int nb_cell, reste;

    nb_lon_cell = m->layer[layer_id]->nlon;
    nb_cell = nb_lon_cell / nb_domain;
    reste = nb_lon_cell % nb_domain;

    if (reste == 0) {
	*from = domain_id * nb_cell;
	*to = *from + nb_cell;
	return;
    }

    if (domain_id < reste) {
	nb_cell++;
	*from = domain_id * nb_cell;
    } else {
	*from = reste * (nb_cell + 1) + (domain_id - reste) * nb_cell;
    }

    *to = *from + nb_cell;
}

void get_domain_bounds(struct mesh_t *m,	
		       int layer_id,		
		       int domain_id,		
		       int nb_domain,		
		       int *from,		
		       int *to)
{
    /** where ends the domain id */
    int nb_lon_cell;
    int nb_cell, reste;
    int virtual_domain_id;

    if (domain_id % 2 == 0) {
    	virtual_domain_id = domain_id;
    } else {
    	virtual_domain_id = fmod(domain_id + nb_domain%2, nb_domain );
    }

    nb_lon_cell = m->layer[layer_id]->nlon;
    nb_cell = nb_lon_cell / nb_domain;
    reste = nb_lon_cell % nb_domain;

    if (reste == 0) {
	*from = virtual_domain_id * nb_cell;
	*to = *from + nb_cell;
	return;
    }

    if (virtual_domain_id < reste) {
	nb_cell++;
	*from = virtual_domain_id * nb_cell;
    } else {
	*from = reste * (nb_cell + 1) + (virtual_domain_id - reste) * nb_cell;
    }

    *to = *from + nb_cell;
}

/** \brief return the number of non empty cell in a given domain.
 *
 @param c the light mesh
 @param m the mesh description
 @param domain the current domain id
 @param nbprocs number of running  process
 **/
int count_nbcells_in_domain(struct cell_info_t ****c,	
			    struct mesh_t *m,		
			    int domain,
			    int nbprocs)
{
    int x, y, z;
    int nb_lat_cell, nb_lon_cell;
    int y_start, y_end;
    int nb = 0;

    for (z = 0; z < m->nlayers; z++) {

	nb_lat_cell = m->layer[z]->nlat;
	nb_lon_cell = m->layer[z]->nlon;

	get_domain_bounds(m, z, domain, nbprocs, &y_start, &y_end);

	for (x = 0; x < nb_lat_cell; x++) {
	    for (y = y_start; y < y_end && y < nb_lon_cell; y++) {
		if (c[z][x][y] && c[z][x][y]->nitems) {
		    nb++;
		}
	    }
	}
    }
    return (nb);
}

/** \brief return the number of non empty cell in a given layer 
 *
 @param c the light mesh
 @param m the mesh description
 @param layer the layer in which to count cells
 **/
int count_nbcells_in_layer(struct cell_info_t ****c,	
			   struct mesh_t *m,
			   int layer)
{
    int x, y, z;
    int nb_lat_cell, nb_lon_cell;
    int nb_cells = 0;

    z = layer;

    nb_lat_cell = m->layer[z]->nlat;
    nb_lon_cell = m->layer[z]->nlon;

    for (x = 0; x < nb_lat_cell; x++) {
	for (y = 0; y < nb_lon_cell; y++) {
	    if (c[z][x][y] && c[z][x][y]->nitems) {
		nb_cells++;
	    }
	}
    }

    return (nb_cells);
}


/*
 * \brief return the number of non empty cell in the light mesh 
 *
 *@param c the light mesh
 *@param m the mesh description
 *
 **/
int count_non_empty_cell_in_cellinfo(struct cell_info_t ****c, struct mesh_t *m)
{
    int nb_lat_cell, nb_lon_cell;
    int x, y, z;
    int nb_cells = 0;

    if (!c)
	return (0);		/* the mesh has none of its cell traversed */
    for (z = 0; z < m->nlayers; z++) {
	nb_lat_cell = m->layer[z]->nlat;
	nb_lon_cell = m->layer[z]->nlon;
	for (x = 0; x < nb_lat_cell; x++) {
	    for (y = 0; y < nb_lon_cell; y++) {
		if (c[z][x][y] && c[z][x][y]->nitems) {
		    nb_cells++;
		}
	    }
	}
    }
    return (nb_cells);
}


/** \brief construct a r2m formated file with cell_info data 
 *
 * Constructs a r2m formated file, which is a full cell information description.
 * Returns number of non empties cells in the light mesh.
 * All non empties cells are printed in the following format, 
 *
 * <B>layer cell_x</B>(lat direction) <B>cell_y</B>(long direction) <BR> 
 * <B>faces hits</B> (6 integers)<BR>
 * <B>blocks hits</B> (list of non null blocks)<BR>
 * <B>nb of rays in the cells</B><BR>
 * <B>rayid latitude longitude depth</B> (first point in the cell) <B>latitude longitude depth</B> (last point in the cell) <B>Length_P Length_S</B><BR>
 * ...<BR>
 * <B>rayid ...</B> (nbrays lines)<BR>
 *
 * Example :
 *\verbatim
  0 10 14
  28 0 0 3 0 0
  17 6 22 39 55 26 42 43 59 38 54 45 62 61 63 47 27 11
  4
  1022 0.092424 -1.445175 2.621287 0.093318 -1.445433 18.345807 16.800429 0.000000
  15378 0.090775 -1.441275 10.000000 0.092157 -1.441702 18.906515 12.801157 0.000000
  30724 0.095867 -1.441776 16.295813 0.095693 -1.441720 18.622799 2.600095 0.000000
  4342 0.094214 -1.441970 19.240258 0.095916 -1.442344 3.566029 25.434061 0.000000
 \endverbatim
 
 @param filename where to write the file
 @param c the light mesh
 @param m the mesh description
 @param domain the current domain id
 @param nbprocs number of running process
 @param rank current process 
 **/
int make_domain_info_file(char *filename,	
			   struct cell_info_t ****c,	
			   struct mesh_t *m,		
			   int domain,			
			   int nbprocs,
			   int rank)
{
    int x, y, z;
    int nb_lat_cell, nb_lon_cell;
    int y_start, y_end;
    int i;
    my_File *fd;
    struct cell_info_t *ci;
    int nb_cell_total;
    int nb_cell_cpt=0;
    int nbwrite;
#ifdef DEBUG
    char *date_stamp;

    date_stamp = get_date_stamp();
    fprintf(stdout, "[%s] make_domain_info_file() starts\n", date_stamp);
    fflush(stdout);
    free(date_stamp);
#endif
    fd = my_fopen(filename, openmode);
    if (!fd) {
	fprintf(stderr, "[process %d] make_domain_info_file() : can't open %s for writting. Exit !\n", 
	        rank, filename);
#ifdef USE_MPI
        MPI_Abort(MPI_COMM_WORLD, 1);
#endif
	exit(1);
    } 

    if (!c) {
	my_fclose(fd);
	return(0);
    }

    my_fprintf(fd, "# format=r2m, generated by %s v%s, domain=%d/%d, mesh=%s\n",
		    PACKAGE, VERSION, domain, nbprocs, m->xml_filename);
    nb_cell_total = count_nbcells_in_domain(c, m, domain, nbprocs);
#ifdef DEBUG
    fprintf(stdout,
	    "*** [process %d]  make_domain_info_file domain[%d] nb_cell_total = %d\n",
	    rank, domain, nb_cell_total);
#endif
    nbwrite = my_fprintf(fd, "%d\n", nb_cell_total);
    if (nbwrite < 0)  {
        fprintf(stderr, "[process %d] make_domain_info_file() : can't write into %s. Exit !\n", 
		rank, filename);
#ifdef USE_MPI
        MPI_Abort(MPI_COMM_WORLD, 1);
#endif
        exit(1);
    }
    

    if (nb_cell_total != 0) {
	for (z = 0; z < m->nlayers; z++) {
	    nb_lat_cell = m->layer[z]->nlat;
	    nb_lon_cell = m->layer[z]->nlon;

	    get_domain_bounds(m, z, domain, nbprocs, &y_start, &y_end);

	    for (x = 0; x < nb_lat_cell; x++) {
		for (y = y_start; y < y_end && y < nb_lon_cell; y++) {
		    ci = c[z][x][y];	/* shorthand */

		    /* allocated cells with nitems=0 are orphan cells : ignored */
		    if (ci == NULL)
			continue;
		    if (ci->nitems == 0)
			continue;

		    nb_cell_cpt++;

		    /* x y z */
		    my_fprintf(fd, "%d %d %d\n", z, x, y);

		    /* face hit */
		    for (i = 0; i < NB_HIT_FACES; i++)
			my_fprintf(fd, "%d ", ci->faces_hit[i]);
		    my_fprintf(fd, "\n");

		    /* block hit */
		    my_fprintf(fd, "%d ", ci->nblocks);
		    for (i = 0; i < ci->nblocks; i++)
			my_fprintf(fd, "%d ", ci->block_hit[i]);
		    my_fprintf(fd, "\n");

		    /* items */
		    my_fprintf(fd, "%d\n", ci->nitems);
		    for (i = 0; i < ci->nitems; i++) {

			my_fprintf(fd, "%ld %f %f %f %f %f %f %f %f\n",
				   ci->item[i].rayid,
				   ci->item[i].in->lat,
				   ci->item[i].in->lon,