ex_utils.c

一个用来实现偏微分方程中网格的计算库

  /* printf("[f] list: %ld, *list: %ld \n", list_ptr, *list_ptr); */
  tlist_ptr = *list_ptr;
  last_head_list_ptr = *list_ptr; /* save last head pointer */
  /* printf("[f] last head list: %ld\n",last_head_list_ptr); */

  while (tlist_ptr )                    /* Walk linked list of file ids/vals */
  {
    if (exoid == tlist_ptr->exo_id )    /* linear search for exodus file id */
    {
      if (tlist_ptr == *list_ptr)       /* Are we at the head of the list? */
        *list_ptr = (*list_ptr)->next;  /*   yes, reset ptr to head of list */
      else                              /*   no, remove this record from chain*/
        last_head_list_ptr->next=tlist_ptr->next;
      free(tlist_ptr);                  /* free up memory */
      break;                            /* Quit if found */
    }
    last_head_list_ptr = tlist_ptr;     /* save last head pointer */
    tlist_ptr = tlist_ptr->next;        /* Loop back if not */
  }

/*  printf("[f] list: %ld *list: %ld **list: %ld\n",
          list_ptr,*list_ptr,(*list_ptr)->value, (*list_ptr)->next); */

}

/*****************************************************************************
*
* ex_get_num_props - get number of properties
*
* entry conditions -
*
* exit conditions -
*
* revision history -
*
*
*****************************************************************************/
int ex_get_num_props (int exoid, ex_entity_type obj_type)
{
  int cntr, varid;
  char var_name[MAX_VAR_NAME_LENGTH+1];
  char  errmsg[MAX_ERR_LENGTH];

  cntr = 0;

  /* loop until there is not a property variable defined; the name of */
  /* the variables begin with an increment of 1 ("xx_prop1") so use cntr+1 */

  while (TRUE)
    {
      switch (obj_type)
	{
	case EX_ELEM_BLOCK:
	  strcpy (var_name, VAR_EB_PROP(cntr+1));
	  break;
	case EX_EDGE_BLOCK:
	  strcpy (var_name, VAR_ED_PROP(cntr+1));
	  break;
	case EX_FACE_BLOCK:
	  strcpy (var_name, VAR_FA_PROP(cntr+1));
	  break;
	case EX_NODE_SET:
	  strcpy (var_name, VAR_NS_PROP(cntr+1));
	  break;
	case EX_EDGE_SET:
	  strcpy (var_name, VAR_ES_PROP(cntr+1));
	  break;
	case EX_FACE_SET:
	  strcpy (var_name, VAR_FS_PROP(cntr+1));
	  break;
	case EX_SIDE_SET:
	  strcpy (var_name, VAR_SS_PROP(cntr+1));
	  break;
	case EX_ELEM_SET:
	  strcpy (var_name, VAR_ELS_PROP(cntr+1));
	  break;
	case EX_ELEM_MAP:
	  strcpy (var_name, VAR_EM_PROP(cntr+1));
	  break;
	case EX_FACE_MAP:
	  strcpy (var_name, VAR_FAM_PROP(cntr+1));
	  break;
	case EX_EDGE_MAP:
	  strcpy (var_name, VAR_EDM_PROP(cntr+1));
	  break;
	case EX_NODE_MAP:
	  strcpy (var_name, VAR_NM_PROP(cntr+1));
	  break;
	default:
	  exerrval = EX_BADPARAM;
	  sprintf(errmsg, "Error: object type %d not supported; file id %d",
		  obj_type, exoid);
	  ex_err("ex_get_prop_names",errmsg,exerrval);
	  return(EX_FATAL);
	}

      if (nc_inq_varid (exoid, var_name, &varid) != NC_NOERR) {
	/*   no variable with this name; return cntr which is now the number of */
	/*   properties for this type of entity */
	return (cntr);
      }
      cntr++;
    }
}

int ex_get_cpu_ws(void)
{
  return(sizeof(float));
}


/* swap - interchange v[i] and v[j] */
static void ex_swap (int v[], int i, int j)
{
  int temp;

  temp = v[i];
  v[i] = v[j];
  v[j] = temp;
}

/*!
 * The following 'indexed qsort' routine is modified from Sedgewicks
 * algorithm It selects the pivot based on the median of the left,
 * right, and center values to try to avoid degenerate cases ocurring
 * when a single value is chosen.  It performs a quicksort on
 * intervals down to the EX_QSORT_CUTOFF size and then performs a final
 * insertion sort on the almost sorted final array.  Based on data in
 * Sedgewick, the EX_QSORT_CUTOFF value should be between 5 and 20.
 *
 * See Sedgewick for further details
 * Define DEBUG_QSORT at the top of this file and recompile to compile
 * in code that verifies that the array is sorted.
 */

#define EX_QSORT_CUTOFF 12

static int ex_int_median3(int v[], int iv[], int left, int right)
{
  int center;
  center = (left + right) / 2;

  if (v[iv[left]] > v[iv[center]])
    ex_swap(iv, left, center);
  if (v[iv[left]] > v[iv[right]])
    ex_swap(iv, left, right);
  if (v[iv[center]] > v[iv[right]])
    ex_swap(iv, center, right);

  ex_swap(iv, center, right-1);
  return iv[right-1];
}

static void ex_int_iqsort(int v[], int iv[], int left, int right)
{
  int pivot;
  int i, j;
  
  if (left + EX_QSORT_CUTOFF <= right) {
    pivot = ex_int_median3(v, iv, left, right);
    i = left;
    j = right - 1;

    for ( ; ; ) {
      while (v[iv[++i]] < v[pivot]);
      while (v[iv[--j]] > v[pivot]);
      if (i < j) {
        ex_swap(iv, i, j);
      } else {
        break;
      }
    }

    ex_swap(iv, i, right-1);
    ex_int_iqsort(v, iv, left, i-1);
    ex_int_iqsort(v, iv, i+1, right);
  }
}

static void ex_int_iisort(int v[], int iv[], int N)
{
  int i,j;
  int ndx = 0;
  int small;
  int tmp;
  
  small = v[iv[0]];
  for (i = 1; i < N; i++) {
    if (v[iv[i]] < small) {
      small = v[iv[i]];
      ndx = i;
    }
  }
  /* Put smallest value in slot 0 */
  ex_swap(iv, 0, ndx);

  for (i=1; i <N; i++) {
    tmp = iv[i];
    for (j=i; v[tmp] < v[iv[j-1]]; j--) {
      iv[j] = iv[j-1];
    }
    iv[j] = tmp;
  }
}

void ex_iqsort(int v[], int iv[], int N)
{
  ex_int_iqsort(v, iv, 0, N-1);
  ex_int_iisort(v, iv, N);

#if defined(DEBUG_QSORT)
  fprintf(stderr, "Checking sort of %d values\n", N+1);
  int i;
  for (i=1; i < N; i++) {
    assert(v[iv[i-1]] <= v[iv[i]]);
  }
#endif
}

/*!
 * Determine whether the new large model storage is being used in this
 * file, or old method. Basically, the difference is whether the
 * coordinates and nodal variables are stored in a blob (xyz
 * components together) or as a variable per component per
 * nodal_variable.
 */
int ex_large_model(int exoid)
{
  if (exoid < 0) {
    /* If exoid not specified, then query is to see if user specified
     * the large model via an environment variable
     */
    char *option = getenv("EXODUS_LARGE_MODEL");
    if (option != NULL) {
      if (option[0] == 'n' || option[0] == 'N') {
        fprintf(stderr,
                "EXODUSII: Small model size selected via EXODUS_LARGE_MODEL environment variable\n");
        return 0;
      } else {
        fprintf(stderr,
                "EXODUSII: Large model size selected via EXODUS_LARGE_MODEL environment variable\n");
        return 1;
      }
    } else {
      return EXODUS_DEFAULT_SIZE; /* Specified in exodusII_int.h */
    }

  } else {
    /* See if the ATT_FILESIZE attribute is defined in the file */
    int file_size = 0; 
    if (nc_get_att_int(exoid, NC_GLOBAL, ATT_FILESIZE, &file_size) != NC_NOERR) {
      /* Variable not found; default is 0 */
      file_size = 0;
    }
    return file_size;
  }
}
  
int ex_get_dimension(int exoid, const char* DIMENSION, const char *label,
                     size_t *count, int *dimid, const char *routine)
{
  char errmsg[MAX_ERR_LENGTH];
  int status;
  
  *count = 0;
  *dimid = -1;
  
  if ((status = nc_inq_dimid(exoid, DIMENSION, dimid)) != NC_NOERR) {
    exerrval = status;
    if (routine != NULL) {
      if (status == NC_EBADDIM) {
        sprintf(errmsg,
                "Warning: no %s defined in file id %d",
                label, exoid);
        ex_err(routine, errmsg,exerrval);
        
      } else {
        sprintf(errmsg,
                "Error: failed to locate number of %s in file id %d",
                label, exoid);
        ex_err(routine,errmsg,exerrval);
      }
    }
    return status;
  }

  if ((status = nc_inq_dimlen (exoid, *dimid, count)) != NC_NOERR) {
    exerrval = status;
    if (routine != NULL) {
      sprintf(errmsg,
              "Error: failed to get number of %s in file id %d",
              label, exoid);
      ex_err(routine,errmsg,exerrval);
      return -1;
    }
  }
  return status;
}

/*! Calculate the number of words of storage required to store the
 * header information.  Total bytes can be obtained by multiplying
 * words by 4.  Size is slightly underestimated since it only
 * considers the bulk data storage...
 */
size_t ex_header_size(int exoid)
{
  const char *routine = NULL;
  int iows = 0;
  size_t ndim = 0;
  size_t num_nodes = 0;
  size_t num_elem = 0;
  size_t num_eblk = 0;
  size_t num_map  = 0;
  size_t num_nset = 0;
  size_t num_sset = 0;
  int mapid;
  int temp;
  
  size_t size = 0;
  /* Get word size (2 = 8-byte reals, 1 = 4-byte reals */
  
  if (nc_flt_code(exoid) == NC_DOUBLE) 
    iows = 2;
  else
    iows = 1;
  
  /* coordinates = (ndim * numnp)*iows + maps  */
  ex_get_dimension(exoid, DIM_NUM_DIM,   "dimension", &ndim,      &temp, routine);
  ex_get_dimension(exoid, DIM_NUM_NODES, "nodes",     &num_nodes, &temp, routine);
  size += iows * ndim * num_nodes;

  /* node maps */
  if (nc_inq_varid(exoid, VAR_NODE_NUM_MAP, &mapid) != -1)
    size += num_nodes;

  ex_get_dimension(exoid, DIM_NUM_NM,   "node maps", &num_map, &temp, routine);
  size += num_map * num_nodes;

  /* Element Data */
  ex_get_dimension(exoid, DIM_NUM_ELEM, "elements",  &num_elem, &temp, routine);
  
  /* Element order map */
  if (nc_inq_varid (exoid, VAR_MAP, &mapid) != -1)
    size += num_elem;
   
  if (nc_inq_varid (exoid, VAR_ELEM_NUM_MAP, &mapid) != -1)
    size += num_elem;

  /* Element map(s) */
  ex_get_dimension(exoid, DIM_NUM_EM,     "element maps",   &num_map, &temp, routine);
  size += num_map * num_elem;

  /* Element Blocks... */
  ex_get_dimension(exoid, DIM_NUM_EL_BLK, "element blocks", &num_eblk, &temp, routine);
  if (num_eblk > 0) {
    /* Allocate storage for element block parameters... */
    int *ids = malloc(num_eblk * sizeof(int));
    int i;

    size += 2*num_eblk; /* status + ids */
    
    ex_get_ids(exoid, EX_ELEM_BLOCK, ids);
    for (i=0; i < num_eblk; i++) {
      int num_elem_this_blk = 0;
      int num_nodes_per_elem = 0;
      int num_attr = 0;
      char elem_type[MAX_STR_LENGTH+1];
      ex_get_elem_block(exoid, ids[i], elem_type, &num_elem_this_blk,
                        &num_nodes_per_elem, &num_attr);
      size += num_elem_this_blk * num_nodes_per_elem;
      size += num_elem_this_blk * num_attr * iows;
    }
    free(ids);
  }
  
  /* Nodesets */
  ex_get_dimension(exoid, DIM_NUM_NS, "nodesets", &num_nset, &temp, routine);
  if (num_nset > 0) {
    /* Allocate storage for nodeset parameters... */
    int *ids = malloc(num_nset * sizeof(int));
    int i;

    size += 2*num_nset; /* Status + ids */
    ex_get_ids(exoid, EX_NODE_SET, ids);
    for (i=0; i < num_nset; i++) {
      int num_nodes_in_set = 0;
      int num_df_in_set = 0;
      ex_get_node_set_param(exoid, ids[i], &num_nodes_in_set, &num_df_in_set);
      size += num_nodes_in_set;
      size += num_df_in_set * iows;
    }
    free(ids);
  }

  /* Sidesets */
  ex_get_dimension(exoid, DIM_NUM_SS, "sidesets", &num_sset, &temp, routine);
  if (num_sset > 0) {
    /* Allocate storage for sideset parameters... */
    int *ids = malloc(num_sset * sizeof(int));
    int i;

    size += 2*num_sset; /* Status + ids */
    ex_get_ids(exoid, EX_SIDE_SET, ids);
    for (i=0; i < num_sset; i++) {
      int num_sides_in_set = 0;
      int num_df_in_set = 0;
      ex_get_side_set_param(exoid, ids[i], &num_sides_in_set, &num_df_in_set);
      size += num_sides_in_set * 2;
      size += num_df_in_set * iows;
    }
    free(ids);
  }

  if (ex_large_model(exoid) == 0 && size > (1<<29)) {

    fprintf(stderr, "ERROR: Size to store header information exceeds 2GB in file id %d\n       File is probably corrupt, rerun with environment variable EXODUS_LARGE_MODEL set.\n", exoid);
  }
  return size;
}