exgssn.c

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

    }
    else if (strncmp(elem_blk_parms[i].elem_type,"SPHERE",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = SPHERE;
      /* set side set node stride */
        elem_blk_parms[i].num_nodes_per_side[0] = 1;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"QUAD",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = QUAD;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 4)
        elem_blk_parms[i].num_nodes_per_side[0] = 2;
      else if (elem_blk_parms[i].num_nodes_per_elem == 5)
        elem_blk_parms[i].num_nodes_per_side[0] = 2;
      else
        elem_blk_parms[i].num_nodes_per_side[0] = 3;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"TRIANGLE",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = TRIANGLE;
      /* set default side set node stride */
      if (ndim == 2)  /* 2d TRIs */
      {
        if (elem_blk_parms[i].num_nodes_per_elem == 3)
          elem_blk_parms[i].num_nodes_per_side[0] = 2;
        else 
          elem_blk_parms[i].num_nodes_per_side[0] = 3;
      }
      else if (ndim == 3)  /* 3d TRIs */
      {
        if (elem_blk_parms[i].num_nodes_per_elem == 3)
          elem_blk_parms[i].num_nodes_per_side[0] = 3;
        else 
          elem_blk_parms[i].num_nodes_per_side[0] = 6;
      }
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"SHELL",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = SHELL;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 2) /* KLUDGE for 2D Shells*/
        elem_blk_parms[i].num_nodes_per_side[0] = 2;
      else if (elem_blk_parms[i].num_nodes_per_elem == 4)
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else
        elem_blk_parms[i].num_nodes_per_side[0] = 8;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"HEX",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = HEX;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 8)  /* 8-node bricks */
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else if (elem_blk_parms[i].num_nodes_per_elem == 9)  /* 9-node bricks */
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else if (elem_blk_parms[i].num_nodes_per_elem == 12)  /* HEXSHELLS */
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else if (elem_blk_parms[i].num_nodes_per_elem == 27)  /* 27-node bricks */
        elem_blk_parms[i].num_nodes_per_side[0] = 9;
      else 
        elem_blk_parms[i].num_nodes_per_side[0] = 8;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"TETRA",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = TETRA;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 4)
        elem_blk_parms[i].num_nodes_per_side[0] = 3;
      else if (elem_blk_parms[i].num_nodes_per_elem == 8)
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else
        elem_blk_parms[i].num_nodes_per_side[0] = 6;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"WEDGE",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = WEDGE;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 6)
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else
        elem_blk_parms[i].num_nodes_per_side[0] = 8;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"PYRAMID",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = PYRAMID;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 5)
        elem_blk_parms[i].num_nodes_per_side[0] = 4;
      else
        elem_blk_parms[i].num_nodes_per_side[0] = 8;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"BEAM",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = BEAM;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 2)
        elem_blk_parms[i].num_nodes_per_side[0] = 2;
      else 
        elem_blk_parms[i].num_nodes_per_side[0] = 3;
    }
    else if ( (strncmp(elem_blk_parms[i].elem_type,"TRUSS",3) == 0) ||
              (strncmp(elem_blk_parms[i].elem_type,"BAR",3) == 0) ||
              (strncmp(elem_blk_parms[i].elem_type,"EDGE",3) == 0) )
    {
      elem_blk_parms[i].elem_type_val = TRUSS;
      /* determine side set node stride */
      if (elem_blk_parms[i].num_nodes_per_elem == 2)
        elem_blk_parms[i].num_nodes_per_side[0] = 2;
      else 
        elem_blk_parms[i].num_nodes_per_side[0] = 3;
    }
    else if (strncmp(elem_blk_parms[i].elem_type,"NULL",3) == 0)
    {
      elem_blk_parms[i].elem_type_val = NULL_ELEMENT;
      elem_blk_parms[i].num_nodes_per_side[0] = 0;
      elem_blk_parms[i].num_elem_in_blk = 0;
    }
    else
    { /* unsupported element type; no problem if no sides specified for
         this element block */
      elem_blk_parms[i].elem_type_val = UNK;
      elem_blk_parms[i].num_nodes_per_side[0] = 0;
    }
    elem_blk_parms[i].elem_blk_id = elem_blk_ids[i];    /* save id */
    elem_ctr += elem_blk_parms[i].num_elem_in_blk;
    elem_blk_parms[i].elem_ctr = elem_ctr;      /* save elem number max */
  }


  /* Allocate space for the ss element to element block parameter index array */
  if (!(ss_parm_ndx=malloc(tot_num_ss_elem*sizeof(int))))
  {
    free(elem_blk_parms);
    free(elem_blk_ids);
    free(ss_elem_ndx);
    free(side_set_side_list);
    free(side_set_elem_list);
    exerrval = EX_MEMFAIL;
    sprintf(errmsg,
"Error: failed to allocate space for side set elem parms index for file id %d"
,
            exoid);
    ex_err("ex_get_side_set_node_list",errmsg,exerrval);
    return (EX_FATAL);
  }


  /* Allocate space for the ss element to node list index array */
  if (!(ss_elem_node_ndx=malloc(tot_num_ss_elem*sizeof(int))))
  {
    free(ss_parm_ndx);
    free(elem_blk_parms);
    free(elem_blk_ids);
    free(ss_elem_ndx);
    free(side_set_side_list);
    free(side_set_elem_list);
    exerrval = EX_MEMFAIL;
    sprintf(errmsg,
"Error: failed to allocate space for side set elem to node index for file id %d",
            exoid);
    ex_err("ex_get_side_set_node_list",errmsg,exerrval);
           
    return (EX_FATAL);
  }

/* Build side set element to node list index and side set element 
   parameter index.
*/
  node_ctr = 0;
  for (i=0;i<tot_num_ss_elem;i++)
  {
    for (j=0; j<num_elem_blks; j++)
    {
      if (elem_blk_parms[j].elem_type_val != NULL_ELEMENT)
      if (side_set_elem_list[i] <= elem_blk_parms[j].elem_ctr)
        break;
    }

    if (j >= num_elem_blks)
    {
      exerrval = EX_BADPARAM;
      sprintf(errmsg,
             "Error: Invalid element number %d found in side set %d in file %d",
              side_set_elem_list[i], side_set_id, exoid);
      free(ss_parm_ndx);
      free(ss_elem_node_ndx);
      free(elem_blk_parms);
      free(elem_blk_ids);
      free(ss_elem_ndx);
      free(side_set_side_list);
      free(side_set_elem_list);
      ex_err("ex_get_side_set_node_list",errmsg,EX_MSG);
      return (EX_FATAL);
    }

    ss_parm_ndx[i] = j; /* assign parameter block index */
    ss_elem_node_ndx[i] = node_ctr;     /* assign node list index */

    /* Update node_ctr (which points to next node in chain */

    /* WEDGEs with 3 node sides (side 4 or 5) are special cases */
    if (elem_blk_parms[j].elem_type_val == WEDGE &&
        (side_set_side_list[i] == 4 || side_set_side_list[i] == 5))
    {
      if (elem_blk_parms[j].num_nodes_per_elem == 6)
        node_ctr += 3;  /* 3 node side */
      else
        node_ctr += 6;  /* 6 node side */
    }
    /* PYRAMIDSs with 3 node sides (sides 1,2,3,4) are also special */
    else if (elem_blk_parms[j].elem_type_val == PYRAMID &&
             (side_set_side_list[i] < 5))
    {
      if (elem_blk_parms[j].num_nodes_per_elem == 5)
        node_ctr += 3;  /* 3 node side */
      else
        node_ctr += 6;  /* 6 node side */
    }
    /* side numbers 3,4,5,6 for SHELLs are also special */
    else if (elem_blk_parms[j].elem_type_val == SHELL &&
        (side_set_side_list[i] > 2 ))
    {
      if (elem_blk_parms[j].num_nodes_per_elem == 4)
        node_ctr += 2;  /* 2 node side */
      else
        node_ctr += 3;  /* 3 node side */
    }
    /* side numbers 3,4,5 for 3d TRIs are also special */
    else if (elem_blk_parms[j].elem_type_val == TRIANGLE &&
             ndim == 3 &&
             side_set_side_list[i] > 2 )
    {
      if (elem_blk_parms[j].num_nodes_per_elem == 3)  /* 3-node TRI */
        node_ctr += 2;  /* 2 node side */
      else   /* 6-node TRI */
        node_ctr += 3;  /* 3 node side */
    }
    else /* all other element types */
      node_ctr += elem_blk_parms[j].num_nodes_per_side[0];
  }

  /* All setup, ready to go ... */

  elem_ctr=0;

  for (j=0; j < tot_num_ss_elem; j++)
  {

    if (side_set_elem_list[ss_elem_ndx[j]] > elem_ctr)
    {
      /* release connectivity array space and get next one */
      if (elem_ctr > 0)
      {
        free(connect);
      }

      /* Allocate space for the connectivity array for new element block */
      if (!(connect=malloc(elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].num_elem_in_blk*
                           elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].num_nodes_per_elem*
                           (int)sizeof(int))))
      {
        free(elem_blk_parms);
        free(elem_blk_ids);
        free(ss_elem_ndx);
        free(ss_elem_node_ndx);
        free(ss_parm_ndx);
        free(side_set_side_list);
        free(side_set_elem_list);
        exerrval = EX_MEMFAIL;
        sprintf(errmsg,
        "Error: failed to allocate space for connectivity array for file id %d",
                exoid);
        ex_err("ex_get_side_set_node_list",errmsg,exerrval);
        return (EX_FATAL);
      }

      /* get connectivity array */
      if (ex_get_elem_conn(
                        exoid,
                        elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].elem_blk_id,
                        connect) == -1)
      {
        free(connect);
        free(elem_blk_parms);
        free(elem_blk_ids);
        free(ss_elem_ndx);
        free(ss_elem_node_ndx);
        free(ss_parm_ndx);
        free(side_set_side_list);
        free(side_set_elem_list);
        sprintf(errmsg,
        "Error: failed to allocate space for connectivity array for file id %d",
                exoid);
        ex_err("ex_get_side_set_node_list",errmsg,exerrval);
        return (EX_FATAL);
      }
      elem_ctr = elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].elem_ctr;
    }
/*  For each side in side set, use the appropriate lookup table to  
      determine the nodes from the connect array. */

    elem_num = side_set_elem_list[ss_elem_ndx[j]]-1;/* element number 0-based*/
    /* calculate the relative element number position in it's block*/

    elem_num_pos = elem_num - 
                  (elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].elem_ctr -
                   elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].num_elem_in_blk);

    /* calculate the beginning of the node list for this element by
         using the ss_elem_node_ndx index into the side_sets_node_index
         and adding the element number position * number of nodes per elem */

    num_nodes_per_elem = 
               elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].num_nodes_per_elem;
    node_pos = ss_elem_node_ndx[ss_elem_ndx[j]];
    connect_offset = num_nodes_per_elem*elem_num_pos;
    side_num = side_set_side_list[ss_elem_ndx[j]]-1;

    switch (elem_blk_parms[ss_parm_ndx[ss_elem_ndx[j]]].elem_type_val)
    {
      case CIRCLE:
      case SPHERE:
      { /* Note: no side-node lookup table is used for this simple case */
        side_set_node_list[node_pos] = connect[connect_offset];
        side_set_node_cnt_list[ss_elem_ndx[j]] = 1;   /* 1 node object */
        break;
      }
      case TRUSS:
      case BEAM:
      { /* Note: no side-node lookup table is used for this simple case */
        side_set_node_list[node_pos] = connect[connect_offset];
        side_set_node_list[node_pos+1] = connect[connect_offset+1];
        side_set_node_cnt_list[ss_elem_ndx[j]] = 2;   /* 2 node object */
        if (num_nodes_per_elem > 2)
        {
          side_set_node_cnt_list[ss_elem_ndx[j]] = 3; /* 3 node object */
          side_set_node_list[node_pos+2] = connect[connect_offset+2];
        }
        break;
      }
      case TRIANGLE:
      {
        if (side_num+1 < 1 || side_num+1 > 5) /* side number range check */
        {
          exerrval = EX_BADPARAM;
          sprintf(errmsg,
                 "Error: Invalid triangle edge number %d in file id %d",
                  side_num+1, exoid);
          ex_err("ex_get_side_set_node_list",errmsg,exerrval);
          free(connect);
          free(elem_blk_parms);
          free(elem_blk_ids);
          free(ss_elem_ndx);
          free(ss_elem_node_ndx);
          free(ss_parm_ndx);
          free(side_set_side_list);
          free(side_set_elem_list);
          return(EX_FATAL);
        }

        if (ndim == 2)   /* 2d TRIs */
        {
          side_set_node_list[node_pos] = 
            connect[connect_offset+tri_table[side_num][0]-1];
          side_set_node_list[node_pos+1] = 
            connect[connect_offset+tri_table[side_num][1]-1];
          side_set_node_cnt_list[ss_elem_ndx[j]] = 2;   /* 2 node object */
          if (num_nodes_per_elem > 3)   /* 6-node TRI  */
          {
            side_set_node_cnt_list[ss_elem_ndx[j]] = 3; /* 3 node object */
            side_set_node_list[node_pos+2] = 
              connect[connect_offset+tri_table[side_num][2]-1];