自学教程：C++ EH函数代码示例

void wr_elem_result_exo(Exo_DB *exo, const char *filename, double ***vector, 		   const int variable_index, const int time_step,		   const double time_value, 		   struct Results_Description *rd){  int error, i;  double local_time_value=time_value;  /* static char *yo = "wr_elem_result_exo"; */  /*   * This file must already exist.   */  exo->cmode = EX_WRITE;  exo->io_wordsize = 0;		/* query */  exo->exoid = ex_open(filename, exo->cmode, &exo->comp_wordsize, 		       &exo->io_wordsize, &exo->version);  EH(exo->exoid, "ex_open");#ifdef DEBUG  fprintf(stderr, "/t/tfilename    = /"%s/"/n", filename);  fprintf(stderr, "/t/tcomp_ws     = %d/n", exo->comp_wordsize);  fprintf(stderr, "/t/tio_wordsize = %d/n", exo->io_wordsize);#endif  error = ex_put_time (exo->exoid, time_step, &local_time_value );  EH(error, "ex_put_time");  /* If the truth table has NOT been set up, this will be really slow... */  for (i = 0; i < exo->num_elem_blocks; i++) {    if (exo->elem_var_tab_exists == TRUE) {      /* Only write out vals if this variable exists for the block */      if (exo->elem_var_tab[i*rd->nev + variable_index] == 1) {	error = ex_put_var(exo->exoid, time_step, EX_ELEM_BLOCK, variable_index+1,				exo->eb_id[i], exo->eb_num_elems[i],				vector[i][variable_index]);	EH(error, "ex_put_var elem");      }    }    else {      /* write it anyway (not really recommended from a performance viewpoint) */      error      = ex_put_var ( exo->exoid,				time_step, EX_ELEM_BLOCK,				variable_index+1, /* Convert to 1 based for						     exodus */				exo->eb_id[i],				exo->eb_num_elems[i],				vector[i][variable_index] );      EH(error, "ex_put_var elem");    }  }  error      = ex_close ( exo->exoid );  EH(error, "ex_close");  return;}

static voidcheck_discontinuous_interp_type(PROBLEM_DESCRIPTION_STRUCT *curr_pd,			        int var_type)         /********************************************************************     *     * check_discontinuous_interp_type     *     * -> Test whether an interpolation specified in the input deck     *    is consistent with a discontinuous interpolation at     *    the interface     ********************************************************************/{  int *v_ptr = curr_pd->v;  int interp_type = curr_pd->i[var_type];  if (v_ptr[var_type] & V_MATSPECIFIC) {    switch (interp_type) {    case I_NOTHING:    case I_Q1:    case I_Q2:    case I_Q2_LSA:    case I_H3:    case I_S2:    case I_B3:    case I_Q3:    case I_Q4:    case I_SP:	fprintf(stderr,"check_interpolation_discontinuous ERROR");	fprintf(stderr," var type to be discontinuous in mat %s/n,",		curr_pd->MaterialName);	fprintf(stderr,"/tbut incompatible interp type specified: %d/n",		interp_type);	EH(-1,"incompatible interp type");	break;	    case I_G0:    case I_G1:    case I_P0:    case I_P1:    case I_Q1_D:    case I_Q2_D:    case I_Q2_D_LSA:    case I_PQ1:    case I_PQ2:         break;    default:     	fprintf(stderr,"check_interpolation_discontinuous ERROR");	fprintf(stderr,"unknown interpolation type/n");	EH(-1,"unknown interpolation type");    }  }}

static voidrearrange_mat_increasing(int var_type, int subvarIndex,			 NODAL_VARS_STRUCT *nv_old,			 NODAL_VARS_STRUCT *nv_new, int *rec_used)          /*****************************************************************       *       * rearrange_mat_increasing():       *       *   Search for the record with the matching variable type and       *   subvar index that  has the lowest MatID field, that is       *   currently unused. It then copies that variable type from       *   the old nodal variable structure to the new.       *****************************************************************/{  int j, mn_min, var_rec, i, num;  VARIABLE_DESCRIPTION_STRUCT *vd;  /*   *  Loop over the number of records in the old structure that   *  have the given variable type.   */  num =  nv_old->Num_Var_Desc_Per_Type[var_type];  if (var_type == MASS_FRACTION) {    num =  nv_old->Num_Var_Desc_Per_Type[var_type] /	   upd->Max_Num_Species_Eqn;  }    for (j = 0; j < num; j++) {     /*     * Search for the record with the matching variable type and     * subvar index that     * has the lowest MatID field, that is currently unused.     */    mn_min = INT_MAX;    var_rec = -1;    for (i = 0; i < nv_old->Num_Var_Desc; i++) {      if (! rec_used[i]) {	vd = nv_old->Var_Desc_List[i];	if ((var_type == vd->Variable_Type) &&	    (subvarIndex == (int) vd->Subvar_Index) &&	    (mn_min > vd->MatID)) {	  var_rec = i;	  mn_min = vd->MatID;	}      }     }    if (var_rec == -1) {      EH(-1,"AUGH! we shouln't be here");    }    rec_used[var_rec] = 1;          /*     * Now fill in fields in the nodal vars struct     */    add_var_to_nv_struct(nv_old->Var_Desc_List[var_rec],			 nv_new);  }}

intgsum_Int(const int value)        /********************************************************************     *     * gsum_Int     *     *   This routine will return the sum of a single integer     *   distributed across the processors.     *     *     *  Return     *  -------     *  Routine returns the sum.     ********************************************************************/{#ifdef PARALLEL  int out_buf, err;  err = MPI_Allreduce((void *) &value, (void *) &out_buf, 1, MPI_INT,		      MPI_SUM, MPI_COMM_WORLD);  if (err != MPI_SUCCESS) {    EH(-1, "gsum_Int: MPI_Allreduce returned an error");  }  return out_buf;#else  return value;#endif}

doublegavg_double(const double value)    /********************************************************************     *     * gavg_double     *     *   This routine will find the average value of an input     *   across all of the processors     *     *     *  Return     *  -------     *  Routine returns the average value     ********************************************************************/{#ifdef PARALLEL  int err;  double out_buf;  err = MPI_Allreduce((void *) &value, (void *) &out_buf, 1, MPI_DOUBLE,	     	      MPI_SUM, MPI_COMM_WORLD);  if (err != MPI_SUCCESS) {    EH(-1, "gavg_double: MPI_Allreduce returned an error");  }  return out_buf / Num_Proc;#else  return value;#endif}

voidretrieve_parameterS(double *lambda, /* PARAMETER VALUE */		  double *x, 	 /* UNKNOWN VECTOR */		  double *xdot,  /* UNKNOWN_DOT VECTOR */		  int sens_type,/*  type of sensitivity variable(BC or MT) */		  int sens_id,		/*  variable id BCID or MT# */		  int sens_flt,		/* data float id or matl prop id */		  int sens_flt2,	/* data float id for UM */		  Comm_Ex *cx,	 /* array of communications structures */		  Exo_DB *exo,	 /* ptr to the finite element mesh database */		  Dpi *dpi)	 /* distributed processing information */{  int ic;  int mn,mpr;  int ibc, idf;  #ifdef DEBUG  static const char yo[]="retrieve_parameterS";  fprintf(stderr, "%s() begins.../n", yo);#endif  if (sens_type == 1) {    ibc = sens_id;    idf = sens_flt;    retrieve_BC_parameter(lambda, ibc, idf, cx, exo, dpi);  	}  else  if (sens_type == 2) {    mn = sens_id;    mpr = sens_flt;      retrieve_MT_parameter(lambda, mn, mpr, cx, exo, dpi);	}  else  if (sens_type == 3) {      ibc = sens_id;     idf = sens_flt;       retrieve_AC_parameter(lambda, ibc, idf, cx, exo, dpi); 	}  else   if (sens_type == 4) {    mn = sens_id;    mpr = sens_flt;    ic = sens_flt2;      retrieve_UM_parameter(lambda, mn, mpr, ic, cx, exo, dpi);  }  else EH(-1, "Bad sens. parameter type!"); }/* END of routine retrieve_parameterS  */

intvariable_description_init(VARIABLE_DESCRIPTION_STRUCT **vd_hdl)    /*******************************************************************     *     * variable_description_init:     *     *  Mallocs and initializes a variable description structure     *  to its default state     *     *  Input     *   *el_hdl => If NULL will allocate a structure. If nonnull it will     *              assume that the structure has already been malloced.     *  Return:     *    Success:           CPC_SUCCESS     *    Interface Failure: CPC_PUB_BAD     *    *el_hdl   => address of the new structure that has just been     *                 malloced.     ******************************************************************/{  if (*vd_hdl == NULL) {    *vd_hdl = smalloc(sizeof(VARIABLE_DESCRIPTION_STRUCT));    if (*vd_hdl == NULL) EH(-1, "variable_description_init");  }  return variable_description_default(*vd_hdl);}

/* * Checksum for SRECORD.  Add all the bytes from the record length field through the data, and take * the ones complement.  This includes the address field, which for SRECORDs can be 2 bytes, 3 bytes or * 4 bytes.  In this case, since the upper bytes of the address will be 0 for records of the smaller sizes, * just add all 4 bytes of the address in all cases. *  * The arguments are: *  *  buf		a pointer to the beginning of the data for the record *  length	the length of the data portion of the record *  chunk_len	the length of the record starting at the address and including the checksum *  chunk_addr	starting address for the data in the record * * Returns an unsigned char with the checksum */static unsigned char srec_csum(unsigned char *buf, unsigned int length, int chunk_len, int chunk_addr) {  int sum = chunk_len + (LO(chunk_addr)) + (HI(chunk_addr)) + (EX(chunk_addr)) + (EH(chunk_addr));  unsigned int i;  for(i = 0; i < length; i++) {    sum += buf[i];  }  return ~sum & 0xff;}

int get_nv_offset_idof(NODAL_VARS_STRUCT *nv, const int varType, 		   const int idof, int subVarType,		   VARIABLE_DESCRIPTION_STRUCT **vd_ptr)    /**************************************************************     *     * get_nv_offset_idof():     *     *      This function returns the offset into the local solution     * vector of the idof'th occurrence of variable type, varType.     * MASS_FRACTION variables types also distinguish between      * subvartypes. It also returns the address of the variable     * description structure pertaining to this unknown.     **************************************************************/{  int i, index, offset = -1, nfound;  int num = nv->Num_Var_Desc_Per_Type[varType];  VARIABLE_DESCRIPTION_STRUCT *vd;#ifdef DEBUG_HKM  if (idof >= num) {    EH(-1, "ERROR");  }#endif  if (varType == MASS_FRACTION) {    for (i = 0, nfound = 0; i < num; i++) {      index = nv->Var_Type_Index[varType][i];      vd = nv->Var_Desc_List[index];      if (vd->Subvar_Index == subVarType) {	if (nfound == idof) {	  if (vd_ptr) *vd_ptr = vd;	  offset = nv->Nodal_Offset[index];	  return offset;	} else {	  nfound++;	}      }    }    EH(-1,"NOT FOUND");  } else {    index = nv->Var_Type_Index[varType][idof];    if (vd_ptr) *vd_ptr = nv->Var_Desc_List[index];        offset = nv->Nodal_Offset[index] + subVarType;  }  return offset;}

voidwr_global_result_exo( Exo_DB *exo, 		      const char *filename, 		      const int time_step, 		      const int ngv, 		      double u[] ){     /*****************************************************************      * write_global_result_exo()       *     -- open/write/close EXODUS II db for all global values      *      * The output EXODUS II database contains the original model      * information with some minor QA and info additions, with new       * global data written.      *       ******************************************************************/  int error;  /*    * This capability is deactivated for parallel processing.   * brkfix doesn't support global variables, when this   * changes this restriction should be removed. TAB 3/2002 */  if( u == NULL ) return ; /* Do nothing if this is NULL */  exo->cmode = EX_WRITE;  exo->io_wordsize = 0;		/* query */  exo->exoid = ex_open(filename, exo->cmode, &exo->comp_wordsize, 		       &exo->io_wordsize, &exo->version);  if (exo->exoid < 0) {    EH(-1,"wr_nodal_result_exo: could not open the output file");  }  error = ex_put_var( exo->exoid, time_step, EX_GLOBAL, 1, 0, ngv, u );  EH(error, "ex_put_var glob_vars");  error = ex_close( exo->exoid);  return;}

voidupdate_user_TP_parameter(double lambda, double *x, double *xdot,                         double *x_AC, Comm_Ex *cx, Exo_DB *exo, Dpi *dpi)/* * This function handles updates to the second (TP) parameter * when LOCA bifurcation tracking algorithms (turning point, * pitchfork, or Hopf are used. It works the same as the above * function, but here lambda is the TP parameter. * * The examples in the above function also apply identically to this * one, and for brevity are not repeated. The standard call to * do_user_update is also the same as above. The main difference is * the variables first_cp and first_tp are set up to indicate which * parameter is being updated and whether this is the first update * call for that parameter. The user needn't worry about this detail. */{    /* Actual function begins here (refer to previous example) */    static int first_tp = 1;    // int first_cp = -1;    // int n = 0;    // int Type, BCID, DFID, MTID, MPID, MDID;    // double value;    /* Declare any additional variables here */    /* If using this function, comment out this line. */    EH(-1, "No user TP continuation conditions entered!");    /* Evaluate any intermediate quantities and/or assign constants here */    /* Enter each continuation condition in sequence in this space */    /* ID's */    /* Value */    /* Update call - copy from example */    /* Done */    first_tp = 0;    return;} /* END of routine update_user_TP_parameter */

void assign_species_desc_suffix(const int species_var_name, char * retn_string)   /***************************************************************************    *    * assign_species_desc_suffix:    *    *  This function assigns a string to the return string, associated with the    *  type of species variable being considered.    *    *  Input    * --------    *  species_var_name: Valid species type    *                   (The valid species types are listed in rf_fem_const.h)    * Output    * --------    *  retn_string:  Has to have an input length of at least 24    ***************************************************************************/{  if (retn_string == NULL) {   EH(-1, "assign_species_desc_suffix: bad interface/n");  }  switch (species_var_name) {  case SPECIES_MASS_FRACTION:      (void) strcpy(retn_string, "Mass Fraction");      break;  case SPECIES_MOLE_FRACTION:      (void) strcpy(retn_string, "Mole Fraction");      break;  case SPECIES_VOL_FRACTION:      (void) strcpy(retn_string, "Volume Fraction");      break;  case SPECIES_DENSITY:      (void) strcpy(retn_string, "Density");      break;  case SPECIES_CONCENTRATION:      (void) strcpy(retn_string, "Concentration");      break;  case SPECIES_CAP_PRESSURE:      (void) strcpy(retn_string, "Capillary Pressure");      break;  case SPECIES_UNDEFINED_FORM:      (void) strcpy(retn_string, "Species Unknown");      break;  default:      (void) strcpy(retn_string, "Mass Fraction(default)");      printf("assign_species_desc_suffix: WARNING unknown form of "	     "species vars: %d/n", species_var_name);  }}

static voiddefine_dimension(const int unit,		 const char *string,		 const int value,		 int *identifier){  int err;  char err_msg[MAX_CHAR_ERR_MSG];  /*   * Dimensions with value zero are problematic, so filter them out.   */  if ( value <= 0 )    {      *identifier = -1;      return;    }#ifdef NETCDF_3    err  = nc_def_dim(unit, string, value, identifier);  if ( err != NC_NOERR )    {      sprintf(err_msg, "nc_def_dim() on %s [<%d] id=%d", string, 		   value, *identifier);      EH(-1, err_msg);    }#endif#ifdef NETCDF_2  err  = ncdimdef(unit, string, value);  sprintf(err_msg, "ncdimdef() on %s [<%d] rtn %d", string, 		 value, err);  EH(err, err_msg);  *identifier = err;#endif    return;}

intdof_lnode_var_type(const int n, const int Element_Type,		   const int proc_node_num, const int var_type,		   PROBLEM_DESCRIPTION_STRUCT *pd_ptr)    /**********************************************************************     *     * dof_lnode_var_type:     *     *  This is a wrapper around dof_lnode_interp_type(). It calculates     *  the number of degrees of freedom located at this local node number,     *  that is actually interpolated on this element (i.e., active on this     *  element)     *  given the interpolation type, given the variable type and the current     *  problem description structure. Within that structure it defines     *  the interpolation to be used for that variable type on that     *  element type.     *  See the dof_lnode_interp_type() description for more information.     *     *  proc_node_num is needed to look up whether this node is on the     *  edge of a domain.     *     *  So, if a variable is located at that node, but is not active for     *  that element, the answer is zero.     *     *  The return value is equal to the number of degrees of freedom at     *  a local node for a variable type corresponding to the current     *  element, given the problem description structure.     *********************************************************************/{  int retn, interp_type, edge;  /*   * Store the interpolation type for the input variable   */  interp_type =  pd_ptr->i[var_type];  /*   * Store whether this node is on an edge of the grid or not   */  edge = (int) Nodes[proc_node_num]->EDGE;  /*   * Now, given the local node number, n, the element type, and the   * interpolation type, return the number of degrees of freedom   */  retn = dof_lnode_interp_type(n, Element_Type, interp_type, edge);  EH(retn, "dof_lnode_var_type: ERROR");    return retn;}

voidReduceBcast_BOR(int *ivec, int length)   /********************************************************************    *    * ReduceBcast_BOR()    *    *  Does a logical bitwize OR operation on a vector of ints    *  distibuted across different processors.    ********************************************************************/{#ifdef PARALLEL  int k, err, *ivec_recv;  if (length <= 0) {    printf(" ReduceBcast_BOR Warning, length = %d/n", length);    return;  }  if (ivec == NULL) {    EH(-1, " ReduceBcast_BOR fatal Interface error");  }  ivec_recv = alloc_int_1(length, INT_NOINIT);  err = MPI_Reduce(ivec, ivec_recv, length, MPI_INT, MPI_BOR, 0,	           MPI_COMM_WORLD);  if (err != MPI_SUCCESS) {    EH(-1, " ReduceBcast_BOR fatal MPI Error");  }  err = MPI_Bcast(ivec_recv, V_LAST, MPI_INT, 0, MPI_COMM_WORLD);  if (err != MPI_SUCCESS) {    EH(-1, " ReduceBcast_BOR fatal MPI Error");  }    for (k = 0; k < V_LAST; k++) {    ivec[k] = ivec_recv[k];  }  safer_free((void **) &ivec_recv);#endif}

void wr_nodal_result_exo(Exo_DB *exo, char *filename, double vector[],		    int variable_index, int time_step, double time_value)     /*****************************************************************      * write_nodal_result_exo()       *     -- open/write/close EXODUS II db for 1 nodal var at one      *        time step.      *      * The output EXODUS II database contains the original model      * information with some minor QA and info additions, with new       * nodal value solution data written.      *       ******************************************************************/{  char err_msg[MAX_CHAR_IN_INPUT];  int error;  exo->cmode = EX_WRITE;  exo->io_wordsize = 0;		/* query */  exo->exoid = ex_open(filename, exo->cmode, &exo->comp_wordsize, 		       &exo->io_wordsize, &exo->version);  if (exo->exoid < 0)    {      sr = sprintf(err_msg, 		   "ex_open() = %d on /"%s/" failure @ step %d, time = %g",		   exo->exoid, filename, time_step, time_value);      EH(-1, err_msg);    }  error      = ex_put_time(exo->exoid, time_step, &time_value);  EH(error, "ex_put_time");  error      = ex_put_var(exo->exoid, time_step, EX_NODAL, variable_index, 1,				exo->num_nodes, vector);  EH(error, "ex_put_var nodal");  error      = ex_close(exo->exoid);  return;}

voidAF_restore_Jacobian_Flag(void)    /*************************************************************************     *     *  AF_restore_Jacobian_Flag     *     *     *************************************************************************/{  if (originalChoice_Jacobian == -23) {    EH(-1,"restore_Jacobian_Action_Flag ERROR: no original choice");  }  af->Assemble_Jacobian = originalChoice_Jacobian;  originalChoice_Jacobian = -23;}

voidget_nv_vd_from_offset(NODAL_VARS_STRUCT *nv, int offset,		      VARIABLE_DESCRIPTION_STRUCT **vd_ptr,  int *idof)    /*     *     */{  int i, sum = 0;  VARIABLE_DESCRIPTION_STRUCT *vd;  for (i = 0; i < nv->Num_Var_Desc; i++) {    vd = nv->Var_Desc_List[i];    if (vd->Ndof + sum > offset) {      *vd_ptr = vd;      *idof = offset - sum;      return;    }    sum += vd->Ndof;   }  EH(-1,"RAN OUT OF vds");}

intnsid2nn ( int psid ){  int nsp;			/* node set pointer for this pointset */  nsp     = match_nsid(psid);	/* node set pointer this node set */  if( nsp == -1 )    {      if( Num_Proc == 1 )       {        EH(nsp, "Failed to match nodeset.");       }      else        return(-1);    }  return( Proc_NS_List[ Proc_NS_Pointers[nsp] ] ); /* global node num */}/* END of routine nsid2nn */

intgminloc_int(const int value, const int local_loc, int *global_minloc)    /********************************************************************     *     * gminloc_int():     *     *   This routine will find the minimum integer value input from     *   all processors. It will also return the location from which     *   the minimum occurred. In case of ties, it will return the     *   location from the lowest numbered processor ID with the     *   min value.     *     *  Output     *  -------     *  *global_minloc = Returns the value of local_loc from the     *                   processor containing the min value.     *     *  Return     *  -------     *  Routine returns the minimum value of "value" input from     *  any of the processors.     ********************************************************************/{#ifdef PARALLEL  int  in_buf[2], out_buf[2], err;  in_buf[0] = value;  in_buf[1] = local_loc;  err = MPI_Allreduce((void *)in_buf, (void *) out_buf, 1, MPI_2INT,		      MPI_MINLOC, MPI_COMM_WORLD);  if (err != MPI_SUCCESS) {    EH(-1, "gminloc_int: MPI_Allreduce returned an error");  }  *global_minloc = out_buf[1];  return out_buf[0];#else  *global_minloc = 0;  return value;#endif}