excn2s.c

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

/*
 * Copyright (c) 2005 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Governement
 * retains certain rights in this software.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 * 
 *     * Redistributions in binary form must reproduce the above
 *       copyright notice, this list of conditions and the following
 *       disclaimer in the documentation and/or other materials provided
 *       with the distribution.  
 * 
 *     * Neither the name of Sandia Corporation nor the names of its
 *       contributors may be used to endorse or promote products derived
 *       from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * 
 */
/*****************************************************************************
*
* excn2s - ex_cvt_nodes_to_sides: convert nodes to sides
*
* entry conditions - 
*   input parameters:
*       int     exoid                   exodus file id
*       int     *num_elem_per_set       number of element per set
*       int     *num_nodes_per_set      number of nodes per set
*       int     *side_sets_elem_index   index array of elements into elem list
*       int     *side_sets_node_index   index array of nodes
*       int     *side_sets_elem_list    array of elements
*       int     *side_sets_node_list    array of nodes
*
* exit conditions - 
*       int     *side_sets_side_list    array of sides/faces
*
* revision history - 
*
*  $Id: excn2s.c 2928 2008-07-11 17:45:07Z friedmud $
*
*****************************************************************************/

#include <ctype.h>
#include <string.h>
#include <stdlib.h>
#include "exodusII.h"
#include "exodusII_int.h"

/*!
 * This routine is designed to take the results from retrieving the ExodusI
 * style concatenated side sets to the Exodus II V 2.0 definition
 * uses the element id to get the  coordinate node list,  element block
 * connectivity, element type to
 * convert the side set node list to a side/face list. 

 *  \param    exoid                   exodus file id
 *  \param    *num_elem_per_set       number of element per set
 *  \param    *num_nodes_per_set      number of nodes per set
 *  \param    *side_sets_elem_index   index array of elements into elem list
 *  \param    *side_sets_node_index   index array of nodes
 *  \param    *side_sets_elem_list    array of elements
 *  \param    *side_sets_node_list    array of nodes
 *  \param[out]    *side_sets_side_list    array of sides/faces

 <b>Algorithm:</b>

\verbatim
  Read elem_block_ids --> elem_blk_id[array]

  Read element block parameters --> elem_blk_parms[array]

  Determine total number of elements in side set by summing num_elem_per_set

  Build side set element to side set node list index --> ss_elem_node_ndx[array]

  For each element in the side_set_elem_list  {
    If Jth element is not in current element block (e.g. J>elem_ctr) {
      get element block parameters (num_elem_in_blk, ...)
      elem_ctr += num_elem_in_blk

      free old connectity array space 
      allocate connectivity array: size=num_elem_in_blk*num_nodes_per_elem
      get connectivity array
    }

    If Jth element is in current element block (e.g. J<=elem_ctr) {
      For each node in element (linear search of up to num_nodes_per_elem) {
        If side set element node[1] == element node[i] {
          Case element type = Hex {
            If side set element node[2] == element node[Hex_table[i,1]]
              Jth side = Hex_table[i,2]

             break
          }
          Case element type = Wedge {
            If side set element node[2] == element node[Wedge_table[i,1]]
              Jth side = Wedge_table[i,2]

            break
          }
        }
      }
    }
  }
\endverbatim

 */

int ex_cvt_nodes_to_sides(int exoid,
                          int *num_elem_per_set,
                          int *num_nodes_per_set,
                          int *side_sets_elem_index,
                          int *side_sets_node_index,
                          int *side_sets_elem_list,
                          int *side_sets_node_list,
                          int *side_sets_side_list)
{
  int i, j, k, m, n;
  int  num_side_sets, num_elem_blks;
  int tot_num_elem = 0, tot_num_ss_elem = 0, elem_num = 0, ndim;
  int *elem_blk_ids, *connect;
  int *ss_elem_ndx, *ss_elem_node_ndx, *ss_parm_ndx;
  int elem_ctr, node_ctr, elem_num_pos;
  int num_elem_in_blk, num_nodes_per_elem, num_node_per_side, num_attr;
  int *same_elem_type, el_type;
  float fdum;
  char *cdum, elem_type[MAX_STR_LENGTH+1];

  struct elem_blk_parm  *elem_blk_parms;

/* node to side translation tables - 
     These tables are used to look up the side number based on the
     first and second node in the side/face list. The side node order
     is found in the original Exodus document, SAND87-2997. The element
     node order is found in the ExodusII document, SAND92-2137. These
     tables were generated by following the right-hand rule for determining
     the outward normal. Note: Only the more complex 3-D shapes require
     these tables, the simple shapes are trivial - the first node found
     is also the side number.
*/

  /*    1     2   3    4                                          node 1 */
  static int shell_table[2][8]  = {
      {2,4, 3,1, 4,2, 1,3},                                    /* node 2 */
      {1,2, 1,2, 1,2, 1,2}                                     /* side # */
  };

  /*    1     2   3    4                                          node 1 */
  static int shell_edge_table[2][8]  = {
      {2,4, 3,1, 4,2, 1,3},                                    /* node 2 */
      {3,6, 4,3, 5,4, 6,5}                                     /* side # */
  };

  /*    1     2   3                                               node 1 */
  static int trishell_table[2][6]  = {
      {2,3, 3,1, 1,2},                                         /* node 2 */
      {1,2, 1,2, 1,2}                                          /* side # */
  };

  /*     1      2      3      4                                   node 1 */
  static int tetra_table[2][12]  = {
      {2,3,4, 1,3,4, 4,1,2, 1,2,3},                            /* node 2 */
      {1,4,3, 4,2,1, 2,3,4, 1,2,3}                             /* side # */
  };

#if 0
  static int wedge_table[2][18]  = {
  /*     1      2      3      4      5      6                     node 1 */
      {2,4,3, 5,1,3, 6,1,2, 1,6,5, 6,2,4, 4,3,5},              /* node 2 */
      {1,3,4, 1,4,2, 2,3,4, 1,3,5, 5,2,1, 5,3,2}               /* side # */
  };
#endif
  
  static int hex_table[2][24]  = {
  /*     1      2      3      4      5      6      7      8       node 1 */
      {4,2,5, 1,3,6, 7,4,2, 3,1,8, 6,8,1, 5,2,7, 8,6,3, 7,5,4},/* node 2 */
      {5,1,4, 5,2,1, 2,3,5, 5,4,3, 6,4,1, 1,2,6, 6,2,3, 3,6,4} /* side # */
  };

  char errmsg[MAX_ERR_LENGTH];

  exerrval = 0; /* clear error code */

  cdum = 0; /* initialize even though it is not used */

/* first check if any side sets are specified */
/* inquire how many side sets have been stored */

  if ((ex_inquire(exoid, EX_INQ_SIDE_SETS, &num_side_sets, &fdum, cdum)) == -1)
  {
    sprintf(errmsg,
           "Error: failed to get number of side sets in file id %d",exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return(EX_FATAL);
  }

  if (num_side_sets == 0)
  {
    sprintf(errmsg,
           "Warning: no side sets defined in file id %d",exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,EX_WARN);
    return(EX_WARN);
  }

  if ((ex_inquire(exoid, EX_INQ_ELEM_BLK, &num_elem_blks, &fdum, cdum)) == -1)
  {
    sprintf(errmsg,
           "Error: failed to get number of element blocks in file id %d",exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return(EX_FATAL);
  }

  if ((ex_inquire(exoid, EX_INQ_ELEM, &tot_num_elem, &fdum, cdum)) == -1)
  {
    sprintf(errmsg,
           "Error: failed to get total number of elements in file id %d",exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return(EX_FATAL);
  }

/* get the dimensionality of the coordinates;  this is necessary to
   distinguish between 2d TRIs and 3d TRIs */

  if ((ex_inquire(exoid, EX_INQ_DIM, &ndim, &fdum, cdum)) == -1)
  {
    sprintf(errmsg,
           "Error: failed to get dimensionality in file id %d",exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return(EX_FATAL);
  }

  /* First count up # of elements in the side sets*/
  for (i=0;i<num_side_sets;i++)
    tot_num_ss_elem += num_elem_per_set[i];

  /* Allocate space for the ss element index array */
  if (!(ss_elem_ndx=malloc(tot_num_ss_elem*sizeof(int))))
  {
    exerrval = EX_MEMFAIL;
    sprintf(errmsg,
 "Error: failed to allocate space for side set elem sort array for file id %d",
            exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return (EX_FATAL);
  }

  /* Sort side set element list into index array  - non-destructive */
  for (i=0;i<tot_num_ss_elem;i++)
    ss_elem_ndx[i] = i; /* init index array to current position */

  ex_iqsort(side_sets_elem_list, ss_elem_ndx,tot_num_ss_elem);


  /* Allocate space for the element block ids */
  if (!(elem_blk_ids=malloc(num_elem_blks*sizeof(int))))
  {
    free(ss_elem_ndx);
    exerrval = EX_MEMFAIL;
    sprintf(errmsg,
        "Error: failed to allocate space for element block ids for file id %d",
            exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return (EX_FATAL);
  }

  if (ex_get_elem_blk_ids(exoid, elem_blk_ids))
  {
    free(elem_blk_ids);
    free(ss_elem_ndx);
    sprintf(errmsg,
           "Error: failed to get element block ids in file id %d",
            exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,EX_MSG);
    return(EX_FATAL);
  } 

  /* Allocate space for the element block params */
  if (!(elem_blk_parms=malloc(num_elem_blks*sizeof(struct elem_blk_parm))))
  {
    free(elem_blk_ids);
    free(ss_elem_ndx);
    exerrval = EX_MEMFAIL;
    sprintf(errmsg,
      "Error: failed to allocate space for element block params for file id %d",
            exoid);
    ex_err("ex_cvt_nodes_to_sides",errmsg,exerrval);
    return (EX_FATAL);
  }
  elem_ctr = 0;
  for (i=0; i<num_elem_blks; i++)
  {
    /* read in an element block parameter */
    if ((ex_get_elem_block (exoid,
                            elem_blk_ids[i],
                            elem_type,
                            &num_elem_in_blk,
                            &num_nodes_per_elem,
                            &num_attr)) == -1)
    {
      free(elem_blk_parms);
      free(elem_blk_ids);
      free(ss_elem_ndx);
      sprintf(errmsg,
             "Error: failed to get element block %d parameters in file id %d",
              elem_blk_ids[i], exoid);
      ex_err("ex_cvt_nodes_to_sides",errmsg,EX_MSG);
      return(EX_FATAL);
    }

    elem_blk_parms[i].num_elem_in_blk = num_elem_in_blk;
    elem_blk_parms[i].num_nodes_per_elem = num_nodes_per_elem;
    elem_blk_parms[i].num_attr = num_attr;