exodusii_io.c

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

    return nodal_var_values;
  }

  // For Writing Solutions

  void ExodusII::create(std::string filename)
  {
    //Store things as doubles
    comp_ws = 8;
    io_ws = 8;
    
    ex_id = exII::ex_create(filename.c_str(), EX_CLOBBER, &comp_ws, &io_ws);
    
    ex_id = exII::ex_open(filename.c_str(),
			  EX_WRITE,
			  &comp_ws,
			  &io_ws,
			  &ex_version);
  
    check_err(ex_id, "Error creating ExodusII mesh file.");
    if (verbose) std::cout << "File created successfully." << std::endl;
  }

  void ExodusII::initialize(std::string title, const MeshBase & mesh)
  {
    num_dim = mesh.mesh_dimension();
    num_nodes = mesh.n_nodes();
    num_elem = mesh.n_elem();
    num_elem_blk = 1;
    num_node_sets = 0;
    num_side_sets = 0;

    ex_err = exII::ex_put_init(ex_id, title.c_str(), num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets);
    
    check_err(ex_err, "Error initializing new Exodus file.");
  }

  void ExodusII::write_nodal_coordinates(const MeshBase & mesh)
  {
    x.resize(num_nodes);
    y.resize(num_nodes);
    z.resize(num_nodes);

    for (/* unsigned */ int i=0; i<num_nodes; ++i)
    {
      x[i]=(*mesh.node_ptr(i))(0);
      y[i]=(*mesh.node_ptr(i))(1);
      z[i]=(*mesh.node_ptr(i))(2);
    }

    ex_err = exII::ex_put_coord(ex_id, &x[0], &y[0], &z[0]);

    check_err(ex_err, "Error writing coordinates to Exodus file.");
  }

  void ExodusII::write_elements(const MeshBase & mesh)
  {
    ExodusII::ElementMaps em;
    const ExodusII::Conversion conv = em.assign_conversion(mesh.elem(0)->type());

    num_nodes_per_elem = mesh.elem(0)->n_nodes();
      
    ex_err = exII::ex_put_elem_block(ex_id, 1, conv.exodus_elem_type().c_str(), num_elem,num_nodes_per_elem, 0);
    check_err(ex_err, "Error writing element block.");

    connect.resize(num_elem*num_nodes_per_elem);

    for (int i=0;i<num_elem;i++)
      {
	Elem * elem = mesh.elem(i);

	for(int j=0; j<num_nodes_per_elem; j++)
	  {
	    const unsigned int connect_index   = (i*num_nodes_per_elem)+j;
	    const unsigned int elem_node_index = conv.get_node_map(j);
	    if (verbose)
	      {
		std::cout << "Exodus node index: " << j
			  << "=LibMesh node index " << elem_node_index << std::endl;
	      }
	    connect[connect_index] = elem->node(elem_node_index)+1;
	  }
      }

    ex_err = exII::ex_put_elem_conn(ex_id, 1, &connect[0]);
    check_err(ex_err, "Error writing element connectivities");
  }

  void ExodusII::initialize_nodal_variables(std::vector<std::string> names)
  {
    num_nodal_vars = names.size();

    ex_err = exII::ex_put_var_param(ex_id, "n", num_nodal_vars);
    check_err(ex_err, "Error setting number of nodal vars.");

    // Dynamically allocate an array of char*.  FIXME: Is there a way
    // to do this without dynamic memory allocation?
    char** var_names = new char*[num_nodal_vars];    

    // For each char* in var_names, allocate enough space and copy from
    // the C++ string into it for passing to the C interface.
    for (int i=0;i<num_nodal_vars;i++)
      {
	var_names[i] = new char[names[i].size()+1];
	std::strcpy(var_names[i], names[i].c_str());
      }

    ex_err = exII::ex_put_var_names(ex_id, "n", num_nodal_vars, var_names);
    check_err(ex_err, "Error setting nodal variable names.");

    // Clean up allocated memory
    for (int i=0;i<num_nodal_vars;i++)
      delete [] var_names[i];
    
    delete [] var_names;
  }

  void ExodusII::write_timestep(int timestep, double time)
  {
    ex_err = exII::ex_put_time(ex_id, timestep, &time);
    check_err(ex_err, "Error writing timestep.");
  }

  void ExodusII::write_nodal_values(int var_id, const std::vector<Number> & values, int timestep)
  {
    ex_err = exII::ex_put_nodal_var(ex_id, timestep, var_id, num_nodes, &values[0]);
    check_err(ex_err, "Error writing nodal values.");
  }  

  // ------------------------------------------------------------
  // ExodusII::Conversion class members
  const ExodusII::Conversion ExodusII::ElementMaps::assign_conversion(const std::string type)
  {
    if ((type == "QUAD4") || (type == "QUAD"))
      return assign_conversion(QUAD4);

    else if (type == "QUAD8")
      return assign_conversion(QUAD8);

    else if (type == "QUAD9")
      return assign_conversion(QUAD9);

    else if ((type == "TRI3") || (type == "TRIANGLE"))
      return assign_conversion(TRI3);

    else if (type == "TRI6")
      return assign_conversion(TRI6);

    else if ((type == "HEX8") || (type == "HEX") || (type=="hex8")) 
      return assign_conversion(HEX8);

    else if (type == "HEX20")
      return assign_conversion(HEX20);

    else if (type == "HEX27")
      return assign_conversion(HEX27);

    else if ((type == "TETRA4") || (type == "TETRA"))
      return assign_conversion(TET4);

    else if (type == "TETRA10")
      return assign_conversion(TET10);

    else if (type == "WEDGE")
      return assign_conversion(PRISM6);

    else if (type == "PYRAMID" || type == "PYRAMID5")
      return assign_conversion(PYRAMID5);

    else
      {
	std::cerr << "ERROR! Unrecognized element type: " << type << std::endl;
	libmesh_error();
      }

    libmesh_error();
  
    const Conversion conv(tri3_node_map, tri_edge_map, TRI3,"TRI3"); // dummy
    return conv;  
  }



  const ExodusII::Conversion ExodusII::ElementMaps::assign_conversion(const ElemType type)
  {
    switch (type)
      {

      case QUAD4:
	{
	  const Conversion conv(quad4_node_map, quad_edge_map, QUAD4, "QUAD4");
	  return conv;
	}

      case QUAD8:
	{
	  const Conversion conv(quad8_node_map, quad_edge_map, QUAD8, "QUAD8");
	  return conv;
	}
      
      case QUAD9:
	{
	  const Conversion conv(quad9_node_map, quad_edge_map, QUAD9, "QUAD9");
	  return conv;
	}
      
      case TRI3:
	{
	  const Conversion conv(tri3_node_map, tri_edge_map, TRI3, "TRI3");
	  return conv;
	}
      
      case TRI6:
	{
	  const Conversion conv(tri6_node_map, tri_edge_map, TRI6, "TRI6");
	  return conv;
	}
      
      case HEX8:
	{
	  const Conversion conv(hex8_node_map, hex_face_map, HEX8, "HEX8");
	  return conv;
	}
      
      case HEX20:
	{
	  const Conversion conv(hex20_node_map, hex_face_map, HEX20, "HEX20");
	  return conv;
	}
      
      case HEX27:
	{
	  const Conversion conv(hex27_node_map, hex27_face_map, HEX27, "HEX27");
	  return conv;
	}
      
      case TET4:
	{
	  const Conversion conv(tet4_node_map, tet_face_map, TET4, "TETRA4");
	  return conv;
	}
      
      case TET10:
	{
	  const Conversion conv(tet10_node_map, tet_face_map, TET10, "TETRA10");
	  return conv;
	}

      case PRISM6:
	{
	  const Conversion conv(prism6_node_map, prism_face_map, PRISM6, "WEDGE");
	  return conv;
	}

      case PYRAMID5:
	{
	  const Conversion conv(pyramid5_node_map, pyramid_face_map, PYRAMID5, "PYRAMID5");
	  return conv;
	}
	
      default:
	libmesh_error();
      }
    
    libmesh_error();
    
    const Conversion conv(tri3_node_map, tri_edge_map, TRI3, "TRI3"); // dummy
    return conv;  
  }
  
    
//} // end anonymous namespace

#endif // #ifdef HAVE_EXODUS_API




// ------------------------------------------------------------
// ExodusII_IO class members
ExodusII_IO::ExodusII_IO (MeshBase& mesh) :
  MeshInput<MeshBase> (mesh),
  MeshOutput<MeshBase> (mesh),
  ex_ptr (NULL),
  _timestep(1),
  _verbose (false)
{
}

ExodusII_IO::~ExodusII_IO ()
{
#ifndef HAVE_EXODUS_API

  std::cerr <<  "ERROR, ExodusII API is not defined.\n"
	    << std::endl;
  libmesh_error();
    
#else

  if(ex_ptr)
    ex_ptr->close();

#endif
}

void ExodusII_IO::read (const std::string& fname)
{
  // This is a serial-only process for now;
  // the Mesh should be read on processor 0 and
  // broadcast later
  libmesh_assert(libMesh::processor_id() == 0);

#ifndef HAVE_EXODUS_API

  std::cerr <<  "ERROR, ExodusII API is not defined.\n"
	    << "Input file " << fname << " cannot be read"
	    << std::endl;
  libmesh_error();
    
#else
  
  // Get a reference to the mesh we are reading
  MeshBase& mesh = MeshInput<MeshBase>::mesh();
  
  // Clear any existing mesh data
  mesh.clear();
  
  libmesh_assert(mesh.mesh_dimension() != 1); // No support for 1D ExodusII meshes
  
#ifdef DEBUG
    this->verbose() = true;
#endif

    // FIXME[JWP]: Why do we create ex_ptr on the heap?
  this->ex_ptr = new ExodusII(this->verbose()); // Instantiate ExodusII interface
  ExodusII & ex = *ex_ptr;
  
  ExodusII::ElementMaps em;     // Instantiate the ElementMaps interface
    
  ex.open(fname.c_str());       // Open the exodus file, if possible
  ex.read_header();             // Get header information from exodus file
  ex.print_header();            // Print header information

  libmesh_assert(static_cast<unsigned int>(ex.get_num_dim()) == mesh.mesh_dimension()); // Be sure number of dimensions
                                                                                // is equal to the number of 
                                                                                // dimensions in the mesh supplied.
  
  ex.read_nodes();                        // Read nodes from the exodus file
  mesh.reserve_nodes(ex.get_num_nodes()); // Reserve space for the nodes.
  
  // Loop over the nodes, create Nodes with local processor_id 0.
  for (int i=0; i<ex.get_num_nodes(); i++)
    mesh.add_point (Point(ex.get_x(i),
			  ex.get_y(i),
			  ex.get_z(i)), i);
  
  libmesh_assert (static_cast<unsigned int>(ex.get_num_nodes()) == mesh.n_nodes());

  ex.read_block_info();                 // Get information about all the blocks
  mesh.reserve_elem(ex.get_num_elem()); // Reserve space for the elements
   

  // Read in the element connectivity for each block.
  int nelem_last_block = 0;

  // Loop over all the blocks
  for (int i=0; i<ex.get_num_elem_blk(); i++)
    {
      // Read the information for block i
      ex.read_elem_in_block (i);
      int subdomain_id = ex.get_block_id(i);

      // Set any relevant node/edge maps for this element
      const std::string type (ex.get_elem_type());
      const ExodusII::Conversion conv = em.assign_conversion(type); 
      
      // Loop over all the faces in this block
      int jmax = nelem_last_block+ex.get_num_elem_this_blk();
      for (int j=nelem_last_block; j<jmax; j++)
	{
	  Elem* elem = Elem::build (conv.get_canonical_type()).release();
	  libmesh_assert (elem);
          elem->subdomain_id() = subdomain_id;
          elem->set_id(j);
	  mesh.add_elem (elem);
	    
	  // Set all the nodes for this element
	  for (int k=0; k<ex.get_num_nodes_per_elem(); k++)
	    {
	      int gi = (j-nelem_last_block)*ex.get_num_nodes_per_elem() + conv.get_node_map(k); // global index 
	      int node_number   = ex.get_connect(gi);             // Global node number (1-based)
	      elem->set_node(k) = mesh.node_ptr((node_number-1)); // Set node number
	                                                          // Subtract 1 since
		                                                  // exodus is internally 1-based
	    }
	}
      
      // running sum of # of elements per block,
      // (should equal total number of elements in the end)
      nelem_last_block += ex.get_num_elem_this_blk();
    }
  libmesh_assert (static_cast<unsigned int>(nelem_last_block) == mesh.n_elem());
  
  // Read in sideset information -- this is useful for applying boundary conditions
  {
    ex.read_sideset_info(); // Get basic information about ALL sidesets
    int offset=0;
    for (int i=0; i<ex.get_num_side_sets(); i++)
      {
	offset += (i > 0 ? ex.get_num_sides_per_set(i-1) : 0); // Compute new offset
	ex.read_sideset (i, offset);
      }
    
    const std::vector<int>& elem_list = ex.get_elem_list();
    const std::vector<int>& side_list = ex.get_side_list();
    const std::vector<int>& id_list   = ex.get_id_list();
    
    for (unsigned int e=0; e<elem_list.size(); e++)
      {
	// Set any relevant node/edge maps for this element
	const ExodusII::Conversion conv =
	  em.assign_conversion(mesh.elem(elem_list[e]-1)->type());
	
	mesh.boundary_info->add_side (elem_list[e]-1,
				      conv.get_side_map(side_list[e]-1),
				      id_list[e]);
      }
  }

//  ex.close();            // Close the exodus file, if possible
#endif
}



void ExodusII_IO::copy_nodal_solution(System& system, std::string nodal_var_name)
{
  #ifndef HAVE_EXODUS_API

  std::cerr <<  "ERROR, ExodusII API is not defined.\n"
	    << std::endl;
  libmesh_error();
    
  #else

  ExodusII & ex = *ex_ptr;

  std::vector<double> time_steps = ex.get_time_steps();

  //For now just read the first timestep (1)
  const std::vector<double> & nodal_values = ex.get_nodal_var_values(nodal_var_name,1);

  //const DofMap & dof_map = system.get_dof_map();

  const unsigned int var_num = system.variable_number(nodal_var_name);

  for (unsigned int i=0; i<nodal_values.size(); ++i)
  {
    const unsigned int dof_index = MeshInput<MeshBase>::mesh().node_ptr(i)->dof_number(system.number(),var_num,0);

    // If the dof_index is local to this processor, set the value
    if ((dof_index >= system.solution->first_local_index()) && (dof_index <  system.solution->last_local_index()))
      system.solution->set (dof_index, nodal_values[i]);
  }

  system.update();
  
  #endif
}

void ExodusII_IO::write_nodal_data (const std::string& fname,
				    const std::vector<Number>& soln,
				    const std::vector<std::string>& names)
{
  #ifndef HAVE_EXODUS_API

  std::cerr <<  "ERROR, ExodusII API is not defined.\n"
	    << std::endl;
  libmesh_error();
    
  #else

  if (libMesh::processor_id() == 0)
  {
    const MeshBase & mesh = MeshOutput<MeshBase>::mesh();

    int num_vars = names.size();
    int num_nodes = mesh.n_nodes();
  
    if(!ex_ptr)
    {
      ex_ptr = new ExodusII(this->verbose());
      ex_ptr->create(fname);
      ex_ptr->initialize(fname,mesh);
      ex_ptr->write_nodal_coordinates(mesh);
      ex_ptr->write_elements(mesh);
      ex_ptr->initialize_nodal_variables(names);
    }

    for(int c=0; c<num_vars; c++)
    {
      std::vector<Number> cur_soln(num_nodes);

      //Copy out this variable's solution
      for(int i=0; i<num_nodes; i++)
	cur_soln[i] = soln[i*num_vars + c];//c*num_nodes+i];
    
      ex_ptr->write_nodal_values(c+1,cur_soln,_timestep);
    }  
  }
  
  #endif
}

void ExodusII_IO::write_timestep (const std::string& fname,
				  const EquationSystems& es,
				  const int timestep,
				  const double time)
{
#ifndef HAVE_EXODUS_API

  std::cerr <<  "ERROR, ExodusII API is not defined.\n"
	    << std::endl;
  libmesh_error();
    
#else

  _timestep=timestep; 
  write_equation_systems(fname,es);

  if (libMesh::processor_id() == 0)
    ex_ptr->write_timestep(timestep, time);

#endif
}