/* * 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,