/* * Ray2mesh : software for geophysicists. * Compute various scores attached to the mesh cells, based on geometric   information that rays bring when they traverse cells. * * 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 <stdlib.h>#include <string.h>#include <unistd.h>#include <popt.h>#include "ray2mesh.h"#include "r2mfile.h"#include "score.h"#include <sparse/sparse.h>#include <mesh/mesh.h>#include <mesh/cellinfo.h>#include <mesh/import.h>#include <mesh/mesh2xml.h>void parse_command_line(int argc,			char **argv,			char **import_arg,			char **mesh_arg, char **output_arg, 			int *do_sco, int *do_sparse,			int *first_layer, int *last_layer){    static char *import_arg_opt;    static char *mesh_arg_opt;    static char *output_arg_opt;    static int begin_layer=-1;    static int end_layer=-1;    int rc;    poptContext mcx;    enum options {	OPT_ZERO=0,			/* must be 0 */	OPT_HELP,	OPT_MESHFILE,	OPT_IMPORT,	OPT_OUTPUT,	OPT_DO_SCO,	OPT_DO_SPARSE,	OPT_BEGIN_LAYER,        OPT_END_LAYER,	OPT_VERSION    };    const struct poptOption options[] = {	{"help", 'h', POPT_ARG_NONE,	 NULL, OPT_HELP,	 "help", "display this help message"},	{"version", 'v', POPT_ARG_NONE,	 NULL, OPT_VERSION,	 "print version number", NULL},	{"meshfile", 'm', POPT_ARG_STRING,	 &mesh_arg_opt, OPT_MESHFILE,	 "mesh file (xml)", "FILE"},	{"input", 'i', POPT_ARG_STRING,	 &import_arg_opt, OPT_IMPORT,	 "xml files to import (coma separated list)", "FILE"},	{"output", 'o', POPT_ARG_STRING,	 &output_arg_opt, OPT_OUTPUT,	 "output basename", "STRING"},	{"sparse", '\0', POPT_ARG_NONE,	 NULL, OPT_DO_SPARSE,	 "produced sparse file", NULL},	{"sco", '\0', POPT_ARG_NONE,	 NULL, OPT_DO_SCO,	 "produced sco file", NULL},	{"begin-layer", 'b', POPT_ARG_INT, 	 &begin_layer, OPT_BEGIN_LAYER,         "Layer number to begin from", "INT (>=0)"},        {"end-layer", 'e', POPT_ARG_INT,         &end_layer, OPT_END_LAYER,         "Layer number to end at", "INT (>=0)"},	{NULL, '\0', 0, NULL, 0}	/* end the list */    };    /* args parsing using popt */    mcx = poptGetContext("merger", argc, (const char **) argv, options, 0);    if (argc == 1) {	poptPrintUsage(mcx, stdout, 0);	exit(1);    }    do {	rc = poptGetNextOpt(mcx);	switch (rc) {	case OPT_HELP:	    poptPrintHelp(mcx, stdout, 0);	    exit(0);	case OPT_MESHFILE:	    *mesh_arg = mesh_arg_opt;	    break;	case OPT_IMPORT:	    *import_arg = import_arg_opt;	    break;	case OPT_OUTPUT:	    *output_arg = output_arg_opt;	    break;	case OPT_DO_SPARSE:	    *do_sparse = 1;	    break;	case OPT_DO_SCO:	    *do_sco = 1;	    break;	case OPT_VERSION:	    printf("merger from %s-%s\n", PACKAGE, VERSION);	    exit(0);	    break;	case POPT_ERROR_BADOPT:	/* error */	    fprintf(stderr, "%s: %s\n", poptBadOption(mcx, 0),		    poptStrerror(rc));	    poptPrintUsage(mcx, stdout, 0);	    exit(1);	}    }    while (rc > 0);    poptFreeContext(mcx);    /* check mesh file */    if (!mesh_arg_opt) {	fprintf(stderr, "missing -m option\n");	exit(1);    } else {	if (access(mesh_arg_opt, R_OK)) {	    perror(mesh_arg_opt);	    exit(1);	}    }    if (!output_arg_opt) {	fprintf(stderr, "basename for output file is not defined !\n");	exit(1);    }    *first_layer = begin_layer;    *last_layer = end_layer;}/** brief create a sparse file from a cellinfo array  **/int merger_write_sparse (char *sparse_file, 			 struct cell_info_t ****c, 			 struct mesh_t *m, 			 long int max_ray_id) {    int i;    int x, y, z;    int nb_lat_cell, nb_lon_cell;  /* number of cells in latitude/longitude */    int nb_cells = 0;    struct cell_info_t *ci;    long int lincellid;    struct sparse_matrix_t *A;    struct coord_z3_t cell_id;        A = new_sparse_matrix( max_ray_id + 1, m->ncells, !SPARSE_COL_LINK);    fprintf(stdout, "new sparse matrix initialized (%ldx%ld) ...",	    max_ray_id + 1, m->ncells);    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) {                    /* allocated cells with nitems=0 are orphan -> ignored */		    ci = c[z][x][y];		    cell_id.x=x; cell_id.y=y; cell_id.z=z;                    nb_cells++;		    /* cellid of the cell in the global mesh */        	    lincellid = linearize_cell_id(&cell_id, m);			    		    for (i = 0; i < c[z][x][y]->nitems; i++) {			sparse_set_value(A, ci->item[i].rayid,					 lincellid, 					 ci->item[i].P_length, NULL);		    }                }            }        }    }    fprintf(stdout, "filled (%.3f%%).\n", 100 * (double)A->nb_item/((double)A->nb_line*(double)A->nb_col));    /* save the sparse matrix to file */    write_sparse_matrix(A, sparse_file);    free_sparse_matrix(A);    return(nb_cells);}/*----------------------------------------------------------------------*//* main (merger utility)                                                *//*----------------------------------------------------------------------*/int main(int argc, char **argv){    int buffsize;    struct cell_info_t ****c; 		/**< the light mesh to aggregrate data into */    struct mesh_t *mesh;      		/**< the corresponding mesh description */    struct mesh_t *imported_mesh = NULL;/**< meshes to import */    char *mesh_file = NULL;    char *importfilename = NULL;    char **xmlfilelist = NULL;    int nb_xmlfile = 0;    char *buff;    char *output_arg;        char *xml_output;    char *r2m_file=NULL;    char *sco_file=NULL;    char *sparse_file=NULL;    int i;    int nb_cell;    int section;    int first_layer, last_layer;    int do_sco=0;    int do_sparse=0;    long int max_ray_id;/* mac osx workaround to compile with lam */ #ifdef USE_MPI    int a =0;    if (a) MPI_Init();#endif    parse_command_line(argc, argv, 		       &importfilename, &mesh_file, &output_arg, 		       &do_sco, &do_sparse,		       &first_layer, &last_layer);    /**************************************/    /* test if we can open file to import */    /**************************************/    if (importfilename) {	xmlfilelist = parse_separated_list(importfilename, ",");	nb_xmlfile = 0;	while (xmlfilelist[nb_xmlfile]) {	    if (access(xmlfilelist[nb_xmlfile], R_OK) == -1) {		perror(xmlfilelist[nb_xmlfile]);		exit(1);	    }	    nb_xmlfile++;	}    } else {	fprintf(stderr, "No file to import ... exiting\n");	exit(0);    }    /****************************/    /* main mesh initialization */    /****************************/    mesh = mesh_init_from_file(mesh_file);    assert(mesh);    c = mesh_cellinfo_alloc(mesh);    assert(c);    /* layer filter */    if (first_layer == -1) { 	first_layer = 0;    }    if (last_layer == -1) { 	last_layer = mesh->nlayers-1;    }    if (first_layer > last_layer) {       fprintf(stderr, "error, check your layer ID !\n");       exit(0);    }    fprintf (stdout, "Working on layers %d to %d\n",  		    first_layer, last_layer);    if (first_layer < 0 || first_layer > mesh->nlayers-1) {       fprintf(stderr, "first layer is out of range [%d,%d]\n", 0, 		       mesh->nlayers-1);       exit(0);    }    if (last_layer < 0 ||last_layer > mesh->nlayers-1) {       fprintf(stderr, "last layer is out of range [%d,%d]\n", 0, 		       mesh->nlayers-1);       exit(0);    }        /***************************/    /* slice importation stuff */    /***************************/    /* clean all xml section in the global mesh */    for (section=0; section < NB_MESH_FILE_FORMAT; section++) {    	mesh_remove_data_entry(mesh, section);    }    for (i = 0; i < nb_xmlfile; i++) {	int j,l;	char data_file[256];        struct mesh_offset_t **offset;        imported_mesh = mesh_init_from_file(xmlfilelist[i]);	offset = compute_mesh_offset(mesh, imported_mesh);	for (l=0; l<mesh->nlayers; l++) {		if (offset[l]->lat !=0 || offset[l]->lon !=0 || offset[l]->z !=0) {			fprintf(stderr, "Incompatibles meshes used %s/%s in layer %d\n", 				mesh_file, xmlfilelist[i], l);			exit(1);		}                if (offset[l]) free(offset[l]);        }	free(offset);	/* copy EVT and RES section in the global mesh */	for (j = 0; j < imported_mesh->data[EVT]->ndatafile; j++) {	    mesh_add_data_filename(mesh, EVT, 			    imported_mesh->data[EVT]->filename[j]);	}	for (j = 0; j < imported_mesh->data[RES]->ndatafile; j++) {	    mesh_add_data_filename(mesh, RES, 			    imported_mesh->data[RES]->filename[j]);	}	/* import data */	for (j = 0; j < imported_mesh->data[R2M]->ndatafile; j++) {	    snprintf(data_file, 256, "%s/%s",		     imported_mesh->data[R2M]->directory,		     imported_mesh->data[R2M]->filename[j]);	    fprintf(stdout, "loading %s into memory ... ", data_file);	    fflush(stdout);	    buff = load_file_to_memory(data_file, &buffsize);	    if (!buff) {		    exit(1);	    }	    nb_cell = import_cell_buffer(c, buff, buffsize, 			    		  first_layer, last_layer, 					  0, 0);	    fprintf(stdout, "size %d, imported %d cells\n", buffsize,		    nb_cell);	    fflush(stdout);	    unload_file_from_memory(buff, buffsize);	}	free(xmlfilelist[i]);	free_mesh(imported_mesh);    }    free(xmlfilelist);    /************/    /* sco file */    /************/    if (do_sco) {       sco_file = (char *) 	    malloc((strlen(output_arg) + strlen(".sco") + 1 ) * sizeof(char));       assert(sco_file);       sprintf(sco_file, "%s.sco", output_arg);       fprintf(stdout, "Writing merged file to %s ... ", sco_file);       fflush(stdout);       compute_score (c, mesh);       nb_cell = mesh_cellinfo_write_sco (sco_file, c, mesh);       fprintf(stdout, "%d cells.\n", nb_cell);    }    /***************/    /* sparse file */    /***************/    if (do_sparse) {       sparse_file = (char *) 	    malloc((strlen(output_arg) + strlen(".sparse") + 1 )*sizeof(char));       assert(sparse_file);       sprintf(sparse_file, "%s.sparse", output_arg);       fprintf(stdout, "Writing merged file to %s ... ", sparse_file);       fflush(stdout);       /* find the max ray id         * FIXE-ME : doesn't work with -b and -e options 	*/       max_ray_id =  mesh_cellinfo_get_max_rayid (c, mesh);       fprintf(stderr, "max_ray_id = %ld\n", max_ray_id);       nb_cell = merger_write_sparse (sparse_file, c, mesh, max_ray_id);       fprintf(stdout, "%d cells.\n", nb_cell);    }    /*****************/    /* r2m file name */    /*****************/    r2m_file = (char *) 	    malloc((strlen(output_arg) + strlen(".r2m") + 1 ) * sizeof(char));    assert(r2m_file);    sprintf(r2m_file, "%s.r2m", output_arg);    /* Write the merged r2m file */    fprintf(stdout, "Writing merged file to %s ... ", r2m_file);    fflush(stdout);    nb_cell = make_domain_info_file(r2m_file, c, mesh, 0, 1, 0);    fprintf(stdout, "%d cells.\n", nb_cell);        /*****************/    /* xml file name */    /*****************/    xml_output = (char *) 	    malloc((strlen(output_arg) + strlen(".xml") + 1 ) * sizeof(char));    assert(xml_output);    sprintf(xml_output, "%s.xml", output_arg);    /* add R2M section in the xml */    if (r2m_file) mesh_add_data_filename(mesh, R2M, r2m_file);    if (sco_file) mesh_add_data_filename(mesh, SCO, sco_file);    if (sparse_file) mesh_add_data_filename(mesh, SPARSE, sparse_file);    /* write new xml */    mesh2xml(mesh, xml_output);    free(xml_output);    free_mesh (mesh);    free(sco_file);    free(r2m_file);    free(sparse_file);    return (0);}