/* * 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 "save_to_disk.h"#include "memusage.h"#include "renumber.h"#define  DEFAULT_SLICE 50			/**< the number of raydata to distribute in one round *//*----------------------------------------------------------------------*//* main (sequential code)                                               *//*----------------------------------------------------------------------*/int main(int argc, char **argv){				/* return the nb of ray computed */    struct raydata_t *raydata;    struct ray_config_t ray_config;    struct ray_filter_t ray_filter;    char *meshfile = NULL;    char *inputdatafile = NULL;    char *out_filename = NULL;		/**< cell_data filename really used (may change with --limit */    char *tmpdir=NULL;			/**< only usefull with mpi version */    char *rayfilter_filename = NULL;    FILE *filter_fd = NULL;        int nb_ray_total = 0;    int nb_ray_computed = 0;    int nb_ray_rejected = 0;        int iterate = 0;			/**< is it a unique ray (0) or an input datafile (1) */    int nbread = 0;    int nberr = 0;        float mem_used_start=0;		/**< memory used at the begin of the raytracing process */    float limit = -1.;    float memsize = 0.;    int slice = 0;			/**< size of ray to read */    int nb_ray_stored = 0;		/**< nb of ray used to dump memory to disk */    int nb_save = 0;	   		/**< nb of times we have dumped memory to disk */    int size;				/**< nb of ray to trace */    char *xml_output;        FILE *sparse_fd;    FILE *res_fd;    FILE *event_fd;    /* catch Ctrl-C signal */    signal(SIGINT, emergency_halt);#ifdef DEBUG    /* remember to set MALLOC_TRACE=file */    mtrace();#endif    /******************/    /* parse cmd line */    /******************/    ray_filter.activated=0;    parse_command_line(argc, argv,		       &ray_config,		       &ray_filter,		       &meshfile, 		       &out_filename, 		       &inputdatafile, 		       &rayfilter_filename, 		       &limit, 		       &output_format,		       &tmpdir);    check_file_access(meshfile);    check_file_access(inputdatafile);    fprintf(stderr, "Raytracing mode is ");       switch (ray_config.iterative_mode) {	case ITERATIVE_MODE_OFF:	    fprintf(stderr, "ITERATIVE_MODE_OFF\n");	    break;	case ITERATIVE_MODE_ANGULAR:	    fprintf(stderr, "ITERATIVE_MODE_ANGULAR\n");	    break;	case ITERATIVE_MODE_ON:	    fprintf(stderr, "ITERATIVE_MODE_ON\n");    }    if (ray_filter.activated) {      fprintf(stderr, "Filter activated, removing rays with :\n");     fprintf(stderr, "\tresidual  > %.2fs\n\tdist epi  > %.2fdeg\n\tnb bundle > %d\n\n",		     ray_filter.residual_max, 		     ray_filter.delta_max, 		     ray_filter.nb_event_min);    } else {	fprintf(stderr, "Filter disable\n");    }    /* open the filter file if needed */    if (rayfilter_filename) {       char *filename;       filename = (char *) malloc(sizeof(char)*(                   strlen(rayfilter_filename) + strlen(".sei") + 1));       assert(filename);       sprintf(filename, "%s.sei", rayfilter_filename);       if ((filter_fd = fopen(filename, "w")) == NULL) {          fprintf(stderr, "Cannot open %s for writing ray filtered.\n",                           filename);	  exit(1);       }       free(filename);    }     /***********************************************/    /* set chunk size value ie. nb of rays to send */    /***********************************************/    if (limit > 0) {       /* if user specified a memory limit */       /* slice size if 1/10 of the *estimated* value */       mem_used_start = get_mem_usage();       slice = limit * 1024 / (10 * ONE_RAY_MEM);       fprintf(stdout, "*** limit memory usage to %.1fM\n", limit);       fprintf(stdout, "*** memory limit specified ... forcing slice to %d rays.\n",			slice);    } else {	slice = DEFAULT_SLICE;	fprintf(stdout, "*** setting slice to %d rays.\n", slice);    }    fflush(stdout);        /* init of the mesh->parameters struct  */    if (!(mesh = mesh_init_from_file(meshfile))) {	exit(1);    }    /* open rays input file */    if (!(fdinput = fopen(inputdatafile, "r"))) {	fprintf(stderr, "Can not open input ray file '%s' ... exiting!\n",			inputdatafile);	exit(1);    }	    fprintf(stdout, "** Using data file '%s' for ray computing **\n\n", 		    inputdatafile);    size = get_number_of_lines(fdinput);        /***************************************/    /* open the sparse, res and event file */    /***************************************/    sparse_fd = NULL;    res_fd = NULL;    event_fd = NULL;    change_to_next_files (mesh, out_filename, 		    	  size, 		          &sparse_fd, &res_fd, &event_fd, 		          0 /*rank*/ , 			  0 /*nb_save*/, 			  1 /*want_sparse_file*/);    /*************************************/    /* iterate over input datafile lines */    /*************************************/    while (1) {	/* read a bunch of lines from input file */	raydata = get_raydata(fdinput, slice, &nbread, &nberr);	nb_ray_total += (nbread + nberr);	fprintf(stdout, "\n*** requesting SLICE=%d / nbread=%d\n",		slice, nbread);	if (!raydata) {	    break;	}	bunch_of_ray(raydata, nbread, 0,	/* Put 0 as offset */		     &ray_config, &ray_filter,		     &nb_ray_computed, &nb_ray_rejected,		     filter_fd, sparse_fd, res_fd, event_fd);	fprintf(stdout,		"*** computed %d/%d rays and rejected %d/%d (%.2f%%)\n",		nb_ray_computed, nb_ray_total, nb_ray_rejected,		nb_ray_total,		(float) nb_ray_rejected / (float) (nb_ray_total) * 100);	/* shows memory usage */	if (limit > 0) {	    memsize = get_mem_usage() - mem_used_start;	    fprintf(stdout, "*** MEM(pid=%d) = %.2fMB\n", getpid(), memsize);	}	iterate = !feof(fdinput);	if (!iterate) {	    fclose(fdinput);	    break;	}	/* save memory to disk */	if (    (limit > 0)              && (memsize > limit)	     && (nb_ray_computed != nb_ray_stored))  	{             save_memory_to_disk(output_format, out_filename, cell_info, mesh, 			     0 /*rank*/, 1 /*nbprocs*/, nb_save);	     nb_ray_stored = nb_ray_computed;	     nb_save++;     	     /* next set of files */	     change_to_next_files (mesh, out_filename, size, 			           &sparse_fd, &res_fd, &event_fd, 				   0 /* rank */, 				   nb_save, 				   1 /*want_sparse_file*/);	     fprintf(stderr, "\n*** memory limit (%f MB) reached after %d rays.\n", 			     limit, nb_ray_computed);	}    }    /*********************************************/    /* save the data remaining in memory to disk */    /*********************************************/    if (    (limit > 0) 	 && (nb_ray_computed != nb_ray_stored))     { 	    save_memory_to_disk(output_format, out_filename, cell_info, mesh, 			         0 /*rank*/, 1 /*nbprocs*/, nb_save);     }         /* close sparse, res and evt files */    fclose(sparse_fd);    fclose(res_fd);    fclose(event_fd);    /* everyone close its filtered ray file (if filtering used) */    if (filter_fd) {	fclose(filter_fd);    }	    fprintf(stdout,	    "*** computed %d/%d rays and rejected %d/%d (%.2f%%)****\n",	    nb_ray_computed, nb_ray_total, nb_ray_rejected, nb_ray_total,	    (float) nb_ray_rejected / (float) (nb_ray_total) * 100);#ifndef RAYTRACING_ONLY    /**************************/    /* output r2m or sco file */    /* all data is in memory  */    /**************************/    if (limit < 0)  {       if (strchr(output_format,'r')) {	    char *r2m_file;	    	    /* r2m */	    r2m_file = construct_filename(out_filename, "r2m", 0, 0);	    fprintf(stdout,		    "*** writing cell data (r2m formated) in %s\n", r2m_file);	    mesh_add_data_filename(mesh, R2M, r2m_file);	    make_domain_info_file(r2m_file, cell_info, mesh, 0, 1, 0);	    free(r2m_file);       } else {	    /* sco */	    char *sco_file;    	    fprintf(stdout, "*** score computing ***\n");    	    compute_score(cell_info, mesh);	    sco_file = construct_filename(out_filename, "sco", 0, 0);	    fprintf(stdout,		    "*** writing cell data (sco formated) in %s\n", sco_file);	    mesh_add_data_filename(mesh, SCO, sco_file);	    mesh_cellinfo_write_sco (sco_file, cell_info, mesh);	    free(sco_file);       }    }    /* renumber sparse file only if all data remains in memory */    /* if not you should use the merger software               */    if (limit < 0) {       char * filename;       filename = mesh->data[SPARSE]->filename[0];       renumber_sparse_file (0, nb_ray_computed, filename);       filename = mesh->data[RES]->filename[0];       renumber_res_file (0, nb_ray_computed, filename);       filename = mesh->data[EVT]->filename[0];       renumber_evt_file (0, nb_ray_computed, filename);    }        /* save the xml enrichied with SCO/R2M/SPARSE/RES sections */    xml_output = (char *) malloc((strlen(out_filename) + strlen(".xml") + 1) *			sizeof(char));    assert(xml_output);    sprintf(xml_output, "%s.xml", out_filename);    mesh2xml(mesh, xml_output);    free(xml_output);#endif        free_velocity_model (ray_config.velocity_model);#ifdef DEBUG    muntrace();#endif    return (0);}