exodusii_io.c

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

       * Maps the Exodus edge numbering for triangles.
       * Useful for reading sideset information.
       */ 
      static const int tri_edge_map[3]; 

      /**
       * Maps the Exodus edge numbering for quadrilaterals.
       * Useful for reading sideset information.
       */
      static const int quad_edge_map[4]; 

      
      
      /**
       * 3D maps.  These define
       * mappings from ExodusII-formatted
       * element numberings.
       */

      /**
       * The Hex8 node map.
       * Use this map for bi-linear
       * hexahedral elements in 3D.
       */
      static const int hex8_node_map[8]; 

      /**
       * The Hex20 node map.
       * Use this map for serendipity
       * hexahedral elements in 3D.
       */
      static const int hex20_node_map[20]; 

      /**
       * The Hex27 node map.
       * Use this map for bi-quadratic
       * hexahedral elements in 3D.
       */
      static const int hex27_node_map[27]; 

      /**
       * The Tet4 node map.
       * Use this map for linear
       * tetrahedral elements in 3D.
       */
      static const int tet4_node_map[4]; 

      /**
       * The Tet10 node map.
       * Use this map for quadratic
       * tetrahedral elements in 3D.
       */
      static const int tet10_node_map[10]; 

      /**
       * The Prism6 node map.
       * Use this map for quadratic
       * tetrahedral elements in 3D.
       */
      static const int prism6_node_map[6]; 

      /**
       * The Pyramid5 node map.
       * Use this map for linear
       * pyramid elements in 3D.
       */
      static const int pyramid5_node_map[5]; 


      /**
       * 3D face maps
       */
      
      /**
       * Maps the Exodus face numbering for general hexahedrals.
       * Useful for reading sideset information.
       */
      static const int hex_face_map[6];
      
      /**
       * Maps the Exodus face numbering for 27-noded hexahedrals.
       * Useful for reading sideset information.
       */
      static const int hex27_face_map[6];
      
      /**
       * Maps the Exodus face numbering for general tetrahedrals.
       * Useful for reading sideset information.
       */
      static const int tet_face_map[4];
      
      /**
       * Maps the Exodus face numbering for general prisms.
       * Useful for reading sideset information.
       */
      static const int prism_face_map[5];
      
      /**
       * Maps the Exodus face numbering for general pyramids.
       * Useful for reading sideset information.
       */
      static const int pyramid_face_map[5];
      

      /**
       * @returns a conversion object given an element type name.
       */
      const Conversion assign_conversion(const std::string type);
      
      /**
       * @returns a conversion object given an element type.
       */
      const Conversion assign_conversion(const ElemType type);
    };


  
  private:

  
    /**
     * All of the \p ExodusII
     * API functions return
     * an \p int error value.
     * This function checks
     * to see if the error has
     * been set, and if it has,
     * prints the error message
     * contained in \p msg.
     */
    void check_err(const int error, const std::string msg);

    /**
     * Prints the message defined
     * in \p msg to \p std::cout.
     * Can be turned off if
     * verbosity is set to 0.
     */
    void message(const std::string msg);

    /**
     * Prints the message defined
     * in \p msg to \p std::cout
     * and appends the number
     * \p i to the end of the
     * message.  Useful for
     * printing messages in loops.
     * Can be turned off if
     * verbosity is set to 0.
     */
    void message(const std::string msg, int i);

    const bool verbose;                  // On/Off message flag
    int   comp_ws;                       // ?
    int   io_ws;                         // ?
    int   ex_id;                         // File identification flag
    int   ex_err;                        // General error flag
    int   num_dim;                       // Number of dimensions in the mesh
    int   num_nodes;                     // Total number of nodes in the mesh
    int   num_elem;                      // Total number of elements in the mesh
    int   num_elem_blk;                  // Total number of element blocks
    int   num_node_sets;                 // Total number of node sets
    int   num_side_sets;                 // Total number of element sets
    int   num_elem_this_blk;             // Number of elements in this block
    int   num_nodes_per_elem;            // Number of nodes in each element
    int   num_attr;                      // Number of attributes for a given block
    int   req_info;                      // Generic required info tag
    int   ret_int;                       // Generic int returned by ex_inquire
    int   num_elem_all_sidesets;         // Total number of elements in all side sets
    std::vector<int> block_ids;          // Vector of the block identification numbers
    std::vector<int> connect;            // Vector of nodes in an element
    std::vector<int> ss_ids;             // Vector of the sideset IDs
    std::vector<int> num_sides_per_set;  // Number of sides (edges/faces) in current set
    std::vector<int> num_df_per_set;     // Number of distribution factors per set
    std::vector<int> elem_list;          // List of element numbers in all sidesets
    std::vector<int> side_list;          // Side (face/edge) number actually on the boundary 
    std::vector<int> id_list;            // Side (face/edge) id number
    float ex_version;                    // Version of Exodus you are using
    float ret_float;                     // Generic float returned by ex_inquire
    std::vector<double> x;               // x locations of node points
    std::vector<double> y;               // y locations of node points
    std::vector<double> z;               // z locations of node points
    char    ret_char;                    // Generic char returned by ex_inquire
    char*   title;                       // Problem title
    char*   elem_type;                   // Type of element in a given block

    //Solution Data
    int num_time_steps;
    std::vector<double> time_steps;
    int num_nodal_vars;
    std::vector<std::string> nodal_var_names;
    std::vector<double> nodal_var_values;
};


  // ------------------------------------------------------------
  // ExodusII::ElementMaps static data

  // 2D node map definitions
  const int ExodusII::ElementMaps::quad4_node_map[4] = {0, 1, 2, 3};
  const int ExodusII::ElementMaps::quad8_node_map[8] = {0, 1, 2, 3, 4, 5, 6, 7};
  const int ExodusII::ElementMaps::quad9_node_map[9] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
  const int ExodusII::ElementMaps::tri3_node_map[3]  = {0, 1, 2};
  const int ExodusII::ElementMaps::tri6_node_map[6]  = {0, 1, 2, 3, 4, 5};

  // 2D edge map definitions
  const int ExodusII::ElementMaps::tri_edge_map[3] = {0, 1, 2};
  const int ExodusII::ElementMaps::quad_edge_map[4] = {0, 1, 2, 3};

  // 3D node map definitions
  const int ExodusII::ElementMaps::hex8_node_map[8]   = {0, 1, 2, 3, 4, 5, 6, 7};
  const int ExodusII::ElementMaps::hex20_node_map[20] = { 0,  1,  2,  3,  4,  5,  6,  7,  8,  9,
							  10, 11, 12, 13, 14, 15, 16, 17, 18, 19};  

// Perhaps an older Hex27 node numbering?  This no longer works.
//const int ExodusII::ElementMaps::hex27_node_map[27] = { 1,  5,  6,  2,  0,  4,  7,  3, 13, 17, 14,  9,  8, 16,
//							  18, 10, 12, 19, 15, 11, 24, 25, 22, 26, 21, 23, 20};

// After trial-and-error, we find a nearly identical mapping with Exodus,
// only 20 and 26 are transposed.
const int ExodusII::ElementMaps::hex27_node_map[27] = {
  // Vertex and mid-edge nodes
  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
  // Mid-face nodes and centroid
  26, 21, 22, 23, 24, 25, 20};
//20  21  22  23  24  25  26 // LibMesh indices


  const int ExodusII::ElementMaps::tet4_node_map[4]   = {0, 1, 2, 3};
  const int ExodusII::ElementMaps::tet10_node_map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

  const int ExodusII::ElementMaps::prism6_node_map[6]   = {0, 1, 2, 3, 4, 5};
  const int ExodusII::ElementMaps::pyramid5_node_map[5] = {0, 1, 2, 3, 4};
  
  // 3D face map definitions
  const int ExodusII::ElementMaps::tet_face_map[4]     = {1, 2, 3, 0};
  const int ExodusII::ElementMaps::hex_face_map[6]     = {1, 2, 3, 4, 0, 5};
  const int ExodusII::ElementMaps::hex27_face_map[6]   = {1, 0, 3, 5, 4, 2};
  const int ExodusII::ElementMaps::prism_face_map[5]   = {-1,-1,-1,-1,-1}; // Not Implemented!
  const int ExodusII::ElementMaps::pyramid_face_map[5] = {-1,-1,-1,-1,-1}; // Not Implemented!


  // ------------------------------------------------------------
  // ExodusII class members
  ExodusII::~ExodusII()
  {
    delete [] title;
    delete [] elem_type;
  }



  void ExodusII::check_err(const int err, const std::string msg)
  {
    if (err < 0)
      {
	std::cout << msg << std::endl;
	libmesh_error();
      }
  }



  void ExodusII::message(const std::string msg)
  {
    if (verbose) std::cout << msg << std::endl;
  }



  void ExodusII::message(const std::string msg, int i)
  {
    if (verbose) std::cout << msg << i << "." << std::endl;
  }



  void ExodusII::open(const char* filename)
  {
    ex_id = exII::ex_open(filename,
			  EX_READ,
			  &comp_ws,
			  &io_ws,
			  &ex_version);
  
    check_err(ex_id, "Error opening ExodusII mesh file.");
    if (verbose) std::cout << "File opened successfully." << std::endl;
  }



  void ExodusII::read_header()
  {
    ex_err = exII::ex_get_init(ex_id,
			       title,
			       &num_dim,
			       &num_nodes,
			       &num_elem,
			       &num_elem_blk,
			       &num_node_sets,
			       &num_side_sets);

    check_err(ex_err, "Error retrieving header info.");

    num_time_steps = inquire(exII::EX_INQ_TIME, "Error retrieving time steps");

    exII::ex_get_var_param(ex_id, "n", &num_nodal_vars);

    message("Exodus header info retrieved successfully.");
  }




  void ExodusII::print_header()
  {
    if (verbose)
      std::cout << "Title: \t" << title << std::endl
		<< "Mesh Dimension: \t"   << num_dim << std::endl
		<< "Number of Nodes: \t" << num_nodes << std::endl
		<< "Number of elements: \t" << num_elem << std::endl
		<< "Number of elt blocks: \t" << num_elem_blk << std::endl
		<< "Number of node sets: \t" << num_node_sets << std::endl
		<< "Number of side sets: \t" << num_side_sets << std::endl;
  }



  void ExodusII::read_nodes()
  {
    x.resize(num_nodes);
    y.resize(num_nodes); 
    z.resize(num_nodes); 

    ex_err = exII::ex_get_coord(ex_id,
				static_cast<void*>(&x[0]),
				static_cast<void*>(&y[0]),
				static_cast<void*>(&z[0]));
  
    check_err(ex_err, "Error retrieving nodal data.");
    message("Nodal data retrieved successfully."); 
  }



  void ExodusII::print_nodes()
  {
    for (int i=0; i<num_nodes; i++)
      {
	std::cout << "(" << x[i] << ", " << y[i] << ", " << z[i] << ")" << std::endl;
      }
  }



  void ExodusII::read_block_info()
  {
    block_ids.resize(num_elem_blk);
    ex_err = exII::ex_get_elem_blk_ids(ex_id,
				       &block_ids[0]); // Get all element block IDs.
    // Usually, there is only one
    // block since there is only
    // one type of element.
    // However, there could be more.

    check_err(ex_err, "Error getting block IDs.");
    message("All block IDs retrieved successfully."); 
  }

  int ExodusII::get_block_id(int block)
  {
    libmesh_assert(static_cast<unsigned int>(block) < block_ids.size());
    
    return block_ids[block];
  }

  void ExodusII::read_elem_in_block(int block)
  {
    ex_err = exII::ex_get_elem_block(ex_id,
				     block_ids[block],
				     elem_type,
				     &num_elem_this_blk,
				     &num_nodes_per_elem,
				     &num_attr);
    check_err(ex_err, "Error getting block info.");
    message("Info retrieved successfully for block: ", block); 
  
  
  
    // Read in the connectivity of the elements of this block
    connect.resize(num_nodes_per_elem*num_elem_this_blk);
    ex_err = exII::ex_get_elem_conn(ex_id,
				    block_ids[block],
				    &connect[0]);
  
    check_err(ex_err, "Error reading block connectivity.");
    message("Connectivity retrieved successfully for block: ", block); 
  }



  void ExodusII::read_sideset_info()
  {
    ss_ids.resize(num_side_sets);
    if (num_side_sets > 0)
      {
	ex_err = exII::ex_get_side_set_ids(ex_id,
					   &ss_ids[0]);
	check_err(ex_err, "Error retrieving sideset information.");
	message("All sideset information retrieved successfully."); 

	// Resize appropriate data structures -- only do this once outside the loop
	num_sides_per_set.resize(num_side_sets);
	num_df_per_set.resize(num_side_sets);
  
	req_info = exII::EX_INQ_SS_ELEM_LEN; // Inquire about the length of the
	// concatenated side sets element list
	ex_err = exII::ex_inquire(ex_id,
				  req_info,
				  &ret_int,
				  &ret_float,
				  &ret_char);
	check_err(ex_err, "Error inquiring about side set element list length.");

	//std::cout << "Value returned by ex_inquire was: " << ret_int << std::endl;
	num_elem_all_sidesets = ret_int;	
	elem_list.resize (num_elem_all_sidesets);
	side_list.resize (num_elem_all_sidesets);
	id_list.resize   (num_elem_all_sidesets);
      }
  }


  void ExodusII::read_sideset(int id, int offset)
  {
    ex_err = exII::ex_get_side_set_param(ex_id,
					 ss_ids[id],
					 &num_sides_per_set[id],
					 &num_df_per_set[id]);
    check_err(ex_err, "Error retrieving sideset parameters.");
    message("Parameters retrieved successfully for sideset: ", id);

    ex_err = exII::ex_get_side_set(ex_id,
				   ss_ids[id],
				   &elem_list[offset],
				   &side_list[offset]);
    check_err(ex_err, "Error retrieving sideset data.");
    message("Data retrieved successfully for sideset: ", id);

    for (int i=0; i<num_sides_per_set[id]; i++)
      id_list[i+offset] = ss_ids[id];
  }



  void ExodusII::print_sideset_info()
  {
    for (int i=0; i<num_elem_all_sidesets; i++)
      {
	std::cout << elem_list[i] << " " << side_list[i] << std::endl;
      }
  }



  void ExodusII::close()
  {
    ex_err = exII::ex_close(ex_id);
    check_err(ex_err, "Error closing Exodus file.");
    message("Exodus file closed successfully."); 
  }

  int ExodusII::inquire(int req_info, std::string error_msg)
  {
    ex_err = exII::ex_inquire(ex_id,
			      req_info,
			      &ret_int,
			      &ret_float,
			      &ret_char);
    
    check_err(ex_err, error_msg);

    return ret_int;
  }

  const std::vector<double>& ExodusII::get_time_steps()
  {
    time_steps.resize(num_time_steps);
    exII::ex_get_all_times(ex_id, &time_steps[0]);
    return time_steps;
  }

  const std::vector<std::string>& ExodusII::get_nodal_var_names()
  {
    //Max of 100 variable names
    char *var_names[100];
    
    nodal_var_names.resize(num_nodal_vars);
    for (int i=0;i<num_nodal_vars;i++)
      var_names[i]=new char[MAX_STR_LENGTH+1];

    exII::ex_get_var_names(ex_id, "n", num_nodal_vars, var_names);

    // Copy the char buffers into strings.  Then delete
    // the dynamically allocated char buffers.
    for (int i=0;i<num_nodal_vars;i++)
      {
	nodal_var_names[i]=var_names[i];
	delete [] var_names[i];
      }

    return nodal_var_names;
  }

  const std::vector<double>& ExodusII::get_nodal_var_values(std::string nodal_var_name, int time_step)
  {
    nodal_var_values.resize(num_nodes);
    
    get_nodal_var_names();

    //See if we can find the variable we are looking for
    unsigned int var_index = 0;
    bool found = false;

    found = nodal_var_names[var_index] == nodal_var_name;
    
    while(!found && var_index < nodal_var_names.size())
    {
      var_index++;
      found = nodal_var_names[var_index] == nodal_var_name;
    }

    if(!found)
    {
      std::cerr << "Unable to locate variable named: " << nodal_var_name << std::endl;
      return nodal_var_values;
    }

    exII::ex_get_nodal_var(ex_id, time_step, var_index+1, num_nodes, &nodal_var_values[0]);