gmv_io.c

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

        else
	  out << (*it)->processor_id()+1 << '\n';
      
      out << '\n';
    }


  // If there are *any* variables at all in the system (including
  // p level, or arbitrary cell-based data)
  // to write, the gmv file needs to contain the word "variable"
  // on a line by itself.
  bool write_variable = false;

  // 1.) p-levels
  if (this->p_levels() && mesh_max_p_level)
    write_variable = true;

  // 2.) solution data
  if ((solution_names != NULL) && (v != NULL))
    write_variable = true;

  // 3.) cell-centered data
  if ( !(this->_cell_centered_data.empty()) )
    write_variable = true;
  
  if (write_variable)
    out << "variable\n";


  // optionally write the p-level information
  if (this->p_levels() && mesh_max_p_level)
    {
      out << "p_level 0\n";

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      for ( ; it != end; ++it)
        if (this->subdivide_second_order())
	  for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
            out << (*it)->p_level() << " ";
        else
          out << (*it)->p_level() << " ";
      out << "\n\n";
    }



  
  // optionally write cell-centered data
  if ( !(this->_cell_centered_data.empty()) )
    {
      std::map<std::string, const std::vector<Real>* >::iterator       it  = this->_cell_centered_data.begin();
      const std::map<std::string, const std::vector<Real>* >::iterator end = this->_cell_centered_data.end();      

      for (; it != end; ++it)
	{
	  // write out the variable name, followed by a zero.
	  out << (*it).first << " 0\n";

	  const std::vector<Real>* the_array = (*it).second;
	  
	  // Loop over active elements, write out cell data.  If second-order cells
	  // are split into sub-elements, the sub-elements inherit their parent's
	  // cell-centered data.
	  MeshBase::const_element_iterator       elem_it  = mesh.active_elements_begin();
	  const MeshBase::const_element_iterator elem_end = mesh.active_elements_end(); 
	  
	  for (; elem_it != elem_end; ++elem_it)
	    {
	      const Elem* e = *elem_it;
	      
	      // Use the element's ID to find the value...
	      libmesh_assert (e->id() < the_array->size());
	      const Real the_value = the_array->operator[](e->id());
	      
	      if (this->subdivide_second_order())
		for (unsigned int se=0; se < e->n_sub_elem(); se++)
		  out << the_value << " ";
	      else
		out << the_value << " ";
	    }
	  
	  out << "\n\n";
	}
    }



  
  // optionally write the data
  if ((solution_names != NULL) &&
      (v != NULL))
    {      
      const unsigned int n_vars = solution_names->size();

      if (!(v->size() == mesh.n_nodes()*n_vars))
	std::cerr << "ERROR: v->size()=" << v->size()
		  << ", mesh.n_nodes()=" << mesh.n_nodes()
		  << ", n_vars=" << n_vars
		  << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
		  << std::endl;
      
      libmesh_assert (v->size() == mesh.n_nodes()*n_vars);

      for (unsigned int c=0; c<n_vars; c++)
	{

#ifdef USE_COMPLEX_NUMBERS

	  // in case of complex data, write _tree_ data sets
	  // for each component

	  // this is the real part
	  out << "r_" << (*solution_names)[c] << " 1\n";
	  
	  for (unsigned int n=0; n<mesh.n_nodes(); n++)
	    out << std::setprecision(10) << (*v)[n*n_vars + c].real() << " ";
	  
	  out << '\n' << '\n';


	  // this is the imaginary part
	  out << "i_" << (*solution_names)[c] << " 1\n";
	  
	  for (unsigned int n=0; n<mesh.n_nodes(); n++)
	    out << std::setprecision(10) << (*v)[n*n_vars + c].imag() << " ";
	  
	  out << '\n' << '\n';

	  // this is the magnitude
	  out << "a_" << (*solution_names)[c] << " 1\n";
	  for (unsigned int n=0; n<mesh.n_nodes(); n++)
	    out << std::setprecision(10)
		<< std::abs((*v)[n*n_vars + c]) << " ";

	  out << '\n' << '\n';

#else

	  out << (*solution_names)[c] << " 1\n";
	  
	  for (unsigned int n=0; n<mesh.n_nodes(); n++)
	    out << std::setprecision(10) << (*v)[n*n_vars + c] << " ";
	  
	  out << '\n' << '\n';

#endif
	}
      
    }
  
  // If we wrote any variables, we have to close the variable section now
  if (write_variable)
    out << "endvars\n";

  
  // end of the file
  out << "\nendgmv\n";
}







void GMVIO::write_binary (const std::string& fname,
                          const std::vector<Number>* vec,
                          const std::vector<std::string>* solution_names)
{
  std::ofstream out (fname.c_str());
  
  libmesh_assert (out.good());

  // get a reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();

  unsigned int mesh_max_p_level = 0;
  
  char buf[80];

  // Begin interfacing with the GMV data file
  {
    // write the nodes
    std::strcpy(buf, "gmvinput");
    out.write(buf, std::strlen(buf));
	      
    std::strcpy(buf, "ieeei4r4");
    out.write(buf, std::strlen(buf));
  }


  
  // write the nodes
  {
    std::strcpy(buf, "nodes   ");
    out.write(buf, std::strlen(buf));

    unsigned int tempint = mesh.n_nodes();
    
    std::memcpy(buf, &tempint, sizeof(unsigned int));

    out.write(buf, sizeof(unsigned int));

    // write the x coordinate
    float *temp = new float[mesh.n_nodes()];
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      temp[v] = static_cast<float>(mesh.point(v)(0));
    out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());

    // write the y coordinate
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      temp[v] = static_cast<float>(mesh.point(v)(1));
    out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());

    // write the z coordinate
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      temp[v] = static_cast<float>(mesh.point(v)(2));
    out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());

    delete [] temp;
  }


  // write the connectivity
  {
    std::strcpy(buf, "cells   ");
    out.write(buf, std::strlen(buf));

    unsigned int tempint = mesh.n_active_elem();
    
    std::memcpy(buf, &tempint, sizeof(unsigned int));
    
    out.write(buf, sizeof(unsigned int));

    MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

    switch (mesh.mesh_dimension())
      {

      case 1:
        for ( ; it != end; ++it)
          {
            mesh_max_p_level = std::max(mesh_max_p_level,
                                        (*it)->p_level());

            for(unsigned se = 0; se < (*it)->n_sub_elem(); ++se)
              {
                std::strcpy(buf, "line    ");
                out.write(buf, std::strlen(buf));
	      
                tempint = 2;
                std::memcpy(buf, &tempint, sizeof(unsigned int));
                out.write(buf, sizeof(unsigned int));

                std::vector<unsigned int> conn;
                (*it)->connectivity(se,TECPLOT,conn);
	      
                out.write(reinterpret_cast<char*>(&conn[0]), sizeof(unsigned int)*tempint);
              }
          }
        break;
	
      case 2:
        for ( ; it != end; ++it)
          {
            mesh_max_p_level = std::max(mesh_max_p_level,
                                        (*it)->p_level());

            for(unsigned se = 0; se < (*it)->n_sub_elem(); ++se)
              {
                std::strcpy(buf, "quad    ");
                out.write(buf, std::strlen(buf));
                tempint = 4;
                std::memcpy(buf, &tempint, sizeof(unsigned int));
                out.write(buf, sizeof(unsigned int));
                std::vector<unsigned int> conn;
                (*it)->connectivity(se,TECPLOT,conn);
                out.write(reinterpret_cast<char*>(&conn[0]), sizeof(unsigned int)*tempint);
              }
          }
        break;
      case 3:
        for ( ; it != end; ++it)
          {
            mesh_max_p_level = std::max(mesh_max_p_level,
                                        (*it)->p_level());

            for(unsigned se = 0; se < (*it)->n_sub_elem(); ++se)
              {
                std::strcpy(buf, "phex8   ");
                out.write(buf, std::strlen(buf));
                tempint = 8;
                std::memcpy(buf, &tempint, sizeof(unsigned int));
                out.write(buf, sizeof(unsigned int));
                std::vector<unsigned int> conn;
                (*it)->connectivity(se,TECPLOT,conn);
                out.write(reinterpret_cast<char*>(&conn[0]), sizeof(unsigned int)*tempint);
              }
          }
        break;
      default:
        libmesh_error();
	
      }
  }
  
  
  
  // optionally write the partition information
  if (this->partitioning())
    {
      std::strcpy(buf, "material");
      out.write(buf, std::strlen(buf));
      
      unsigned int tmpint = mesh.n_processors();
      std::memcpy(buf, &tmpint, sizeof(unsigned int));
      out.write(buf, sizeof(unsigned int));

      tmpint = 0; // IDs are cell based
      std::memcpy(buf, &tmpint, sizeof(unsigned int));
      out.write(buf, sizeof(unsigned int));


      for (unsigned int proc=0; proc<mesh.n_processors(); proc++)
        {
          std::sprintf(buf, "proc_%d", proc);
          out.write(buf, 8);
        }

      std::vector<unsigned int> proc_id (mesh.n_active_elem());
      
      unsigned int n=0;
      
      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      for ( ; it != end; ++it)
        for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
          proc_id[n++] = (*it)->processor_id()+1;
      
      
      out.write(reinterpret_cast<char *>(&proc_id[0]),
                sizeof(unsigned int)*proc_id.size());
    }

  // If there are *any* variables at all in the system (including
  // p level, or arbitrary cell-based data)
  // to write, the gmv file needs to contain the word "variable"
  // on a line by itself.
  bool write_variable = false;

  // 1.) p-levels
  if (this->p_levels() && mesh_max_p_level)
    write_variable = true;

  // 2.) solution data
  if ((solution_names != NULL) && (vec != NULL))
    write_variable = true;

  //   // 3.) cell-centered data - unsupported
  //   if ( !(this->_cell_centered_data.empty()) )
  //     write_variable = true;

  if (write_variable)
    {
      std::strcpy(buf, "variable");
      out.write(buf, std::strlen(buf));
    }

  // optionally write the partition information
  if (this->p_levels() && mesh_max_p_level)
    {
      unsigned int n_floats = mesh.n_active_elem();
      for (unsigned int i=0; i != mesh.mesh_dimension(); ++i)
        n_floats *= 2;

      float *temp = new float[n_floats];

      std::strcpy(buf, "p_level");
      out.write(buf, std::strlen(buf));

      unsigned int tempint = 0; // p levels are cell data

      std::memcpy(buf, &tempint, sizeof(unsigned int));
      out.write(buf, sizeof(unsigned int));

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
      unsigned int n=0;

      for (; it != end; ++it)
        for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
          temp[n++] = static_cast<float>( (*it)->p_level() );

      out.write(reinterpret_cast<char *>(temp),
                sizeof(float)*n_floats);

      delete [] temp;
    }

  
   // optionally write cell-centered data
   if ( !(this->_cell_centered_data.empty()) )
     {
       std::cerr << "Cell-centered data not (yet) supported in binary I/O mode!" << std::endl;
       
//        std::map<std::string, const std::vector<Real>* >::iterator       it  = this->_cell_centered_data.begin();
//        const std::map<std::string, const std::vector<Real>* >::iterator end = this->_cell_centered_data.end();      

//        for (; it != end; ++it)
//  	{
//  	  // Write out the variable name ...
//  	  std::strcpy(buf, (*it).first.c_str());
//  	  out.write(buf, std::strlen(buf));
	  
//  	  // ... followed by a zero.
//  	  unsigned int tempint = 0; // 0 signifies cell data
//  	  std::memcpy(buf, &tempint, sizeof(unsigned int));
//  	  out.write(buf, sizeof(unsigned int));

//  	  // Get a pointer to the array of cell-centered data values
//  	  const std::vector<Real>* the_array = (*it).second;

// 	  // Since the_array might contain zeros (for inactive elements) we need to
// 	  // make a copy of it containing just values for active elements.
// 	  const unsigned int n_floats = mesh.n_active_elem() * (1<<mesh.mesh_dimension());
// 	  float *temp = new float[n_floats];

// 	  MeshBase::const_element_iterator       elem_it  = mesh.active_elements_begin();
// 	  const MeshBase::const_element_iterator elem_end = mesh.active_elements_end(); 
// 	  unsigned int n=0;

// 	  for (; elem_it != elem_end; ++elem_it)
// 	    {
// 	      // If there's a seg-fault, it will probably be here!
// 	      const float the_value = static_cast<float>(the_array->operator[]((*elem_it)->id()));
	      
// 	      for (unsigned int se=0; se<(*elem_it)->n_sub_elem(); se++)
// 		temp[n++] = the_value;
// 	    }

	  
//  	  // Write "the_array" directly to the file
//  	  out.write(reinterpret_cast<char *>(temp),
//  		    sizeof(float)*n_floats);

// 	  delete [] temp;
//  	}
     }

  
  
  
  // optionally write the data
  if ((solution_names != NULL) &&
      (vec != NULL))
    {
      float *temp = new float[mesh.n_nodes()];

      const unsigned int n_vars = solution_names->size();
      
      for (unsigned int c=0; c<n_vars; c++)
        {