gmv_io.c

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

#ifdef USE_COMPLEX_NUMBERS
          // for complex data, write three datasets


          // Real part
          std::strcpy(buf, "r_");
          out.write(buf, 2);
          std::strcpy(buf, (*solution_names)[c].c_str());
          out.write(buf, 6);
	  
          unsigned int tempint = 1; // always do nodal data
          std::memcpy(buf, &tempint, sizeof(unsigned int));
          out.write(buf, sizeof(unsigned int));
	  
          for (unsigned int n=0; n<mesh.n_nodes(); n++)
            temp[n] = static_cast<float>( (*vec)[n*n_vars + c].real() );
	  
          out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());


          // imaginary part
          std::strcpy(buf, "i_");
          out.write(buf, 2);
          std::strcpy(buf, (*solution_names)[c].c_str());
          out.write(buf, 6);
	  
          std::memcpy(buf, &tempint, sizeof(unsigned int));
          out.write(buf, sizeof(unsigned int));
	  
          for (unsigned int n=0; n<mesh.n_nodes(); n++)
            temp[n] = static_cast<float>( (*vec)[n*n_vars + c].imag() );
	  
          out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());

          // magnitude
          std::strcpy(buf, "a_");
          out.write(buf, 2);
          std::strcpy(buf, (*solution_names)[c].c_str());
          out.write(buf, 6);
	  
          std::memcpy(buf, &tempint, sizeof(unsigned int));
          out.write(buf, sizeof(unsigned int));
	  
          for (unsigned int n=0; n<mesh.n_nodes(); n++)
            temp[n] = static_cast<float>(std::abs((*vec)[n*n_vars + c]));
	  
          out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());

#else


          std::strcpy(buf, (*solution_names)[c].c_str());
          out.write(buf, 8);
	  
          unsigned int tempint = 1; // always do nodal data
          std::memcpy(buf, &tempint, sizeof(unsigned int));
          out.write(buf, sizeof(unsigned int));
	  
          for (unsigned int n=0; n<mesh.n_nodes(); n++)
            temp[n] = static_cast<float>((*vec)[n*n_vars + c]);
	  
          out.write(reinterpret_cast<char *>(temp), sizeof(float)*mesh.n_nodes());


#endif

	  
        }
    
      delete [] temp;
      
    }

  // If we wrote any variables, we have to close the variable section now
  if (write_variable)
    {
      std::strcpy(buf, "endvars ");
      out.write(buf, std::strlen(buf));
    }

  // end the file
  std::strcpy(buf, "endgmv  ");
  out.write(buf, std::strlen(buf));
}









void GMVIO::write_discontinuous_gmv (const std::string& name,
				     const EquationSystems& es,
				     const bool write_partitioning) const
{
  std::vector<std::string> solution_names;
  std::vector<Number>      v;

  // Get a reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();
  
  es.build_variable_names  (solution_names);
  es.build_discontinuous_solution_vector (v);
  
  if (mesh.processor_id() != 0) return;
  
  std::ofstream out(name.c_str());
  
  libmesh_assert (out.good());

  // Begin interfacing with the GMV data file
  {

    // write the nodes    
    out << "gmvinput ascii" << std::endl << std::endl;
    
    // Compute the total weight
    {
      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      unsigned int tw=0;
      
      for ( ; it != end; ++it)
	tw += (*it)->n_nodes();
      
      out << "nodes " << tw << std::endl;
    }
    


    // Write all the x values
    {
      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 n=0; n<(*it)->n_nodes(); n++)
	  out << (*it)->point(n)(0) << " ";

      out << std::endl;
    }
    
    
    // Write all the y values
    {
      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 n=0; n<(*it)->n_nodes(); n++)
	  out << (*it)->point(n)(1) << " ";

      out << std::endl;
    }
    
    
    // Write all the z values
    {
      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 n=0; n<(*it)->n_nodes(); n++)
	  out << (*it)->point(n)(2) << " ";

      out << std::endl << std::endl;
    }
  }


  
  {
    // write the connectivity
    
    out << "cells " << mesh.n_active_elem() << std::endl;

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

    unsigned int nn=1;
    
    switch (mesh.mesh_dimension())
      {
      case 1:
	{
	  for ( ; it != end; ++it)
	    for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
	      {
		if (((*it)->type() == EDGE2) ||
		    ((*it)->type() == EDGE3) ||
		    ((*it)->type() == EDGE4)
#ifdef ENABLE_INFINITE_ELEMENTS
		    || ((*it)->type() == INFEDGE2)
#endif
		    )
		  {
		    out << "line 2" << std::endl;
		    for (unsigned int i=0; i<(*it)->n_nodes(); i++)
		      out << nn++ << " ";
		    
		  }
		else
		  {
		    libmesh_error();
		  }
		
		out << std::endl;
	      }
	  
	  break;
	}
	
      case 2:
	{
	  for ( ; it != end; ++it)
	    for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
	      {
		if (((*it)->type() == QUAD4) ||
		    ((*it)->type() == QUAD8) || // Note: QUAD8 will be output as one central quad and
			                        // four surrounding triangles (though they will be written
			                        // to GMV as QUAD4s).  
		    ((*it)->type() == QUAD9)
#ifdef ENABLE_INFINITE_ELEMENTS
		    || ((*it)->type() == INFQUAD4)
		    || ((*it)->type() == INFQUAD6)
#endif
		    )
		  {
		    out << "quad 4" << std::endl;
		    for (unsigned int i=0; i<(*it)->n_nodes(); i++)
		      out << nn++ << " ";
		    
		  }
                else if (((*it)->type() == TRI3) ||
                         ((*it)->type() == TRI6))
                  {
                    out << "tri 3" << std::endl;
		    for (unsigned int i=0; i<(*it)->n_nodes(); i++)
		      out << nn++ << " ";
		    
                  }
		else
		  {
		    libmesh_error();
		  }
		
		out << std::endl;
	      }
	  
	  break;
	}
	
	
      case 3:
	{
	  for ( ; it != end; ++it)
	    for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
	      {
		if (((*it)->type() == HEX8) ||
		    ((*it)->type() == HEX20) ||
		    ((*it)->type() == HEX27)
#ifdef ENABLE_INFINITE_ELEMENTS
		    || ((*it)->type() == INFHEX8)
		    || ((*it)->type() == INFHEX16)
		    || ((*it)->type() == INFHEX18)
#endif
		    )
		  {
		    out << "phex8 8" << std::endl;
		    for (unsigned int i=0; i<(*it)->n_nodes(); i++)
		      out << nn++ << " ";
		  }
                else if (((*it)->type() == PRISM6) ||
                         ((*it)->type() == PRISM15) ||
                         ((*it)->type() == PRISM18)
#ifdef ENABLE_INFINITE_ELEMENTS
		         || ((*it)->type() == INFPRISM6)
		         || ((*it)->type() == INFPRISM12)
#endif
			 )
                  {
                    out << "pprism6 6" << std::endl;
		    for (unsigned int i=0; i<(*it)->n_nodes(); i++)
		      out << nn++ << " ";
                  }
                else if (((*it)->type() == TET4) ||
                         ((*it)->type() == TET10))
                  {
                    out << "tet 4" << std::endl;
		    for (unsigned int i=0; i<(*it)->n_nodes(); i++)
		      out << nn++ << " ";
                  }
		else
		  {
		    libmesh_error();
		  }
		
		out << std::endl;
	      }
	  
	  break;
	}
	
      default:
	libmesh_error();
      }
    
    out << std::endl;
  }
  

  
  // optionally write the partition information
  if (write_partitioning)
    {
      out << "material "
	  << mesh.n_processors()
	  << " 0"<< std::endl;

      for (unsigned int proc=0; proc<mesh.n_processors(); proc++)
	out << "proc_" << proc << std::endl;
      
      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      for ( ; it != end; ++it)
	out << (*it)->processor_id()+1 << std::endl;
      
      out << std::endl;
    }


  // Writing cell-centered data is not yet supported in discontinuous GMV files.
  if ( !(this->_cell_centered_data.empty()) )
    {
      std::cerr << "Cell-centered data not (yet) supported for discontinuous GMV files!" << std::endl;
    }

  
  
  // write the data
  {
    const unsigned int n_vars = solution_names.size();
      
    //    libmesh_assert (v.size() == tw*n_vars);

    out << "variable" << std::endl;


    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" << std::endl;
        {
          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 n=0; n<(*it)->n_nodes(); n++)
              out << std::setprecision(10) << v[(n++)*n_vars + c].real() << " ";	    
        }	  	  
        out << std::endl << std::endl;


        // this is the imaginary part
        out << "i_" << solution_names[c] << " 1" << std::endl;	  
        {
          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 n=0; n<(*it)->n_nodes(); n++)
              out << std::setprecision(10) << v[(n++)*n_vars + c].imag() << " ";	    
        }	  
        out << std::endl << std::endl;

        // this is the magnitude
        out << "a_" << solution_names[c] << " 1" << std::endl;
        {
          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 n=0; n<(*it)->n_nodes(); n++)
              out << std::setprecision(10)
                  << std::abs(v[(n++)*n_vars + c]) << " ";	    
        }
        out << std::endl << std::endl;

#else

        out << solution_names[c] << " 1" << std::endl;
        {
          MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
          const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
	    
          unsigned int nn=0;
	  
          for ( ; it != end; ++it)
            for (unsigned int n=0; n<(*it)->n_nodes(); n++)
              out << std::setprecision(10) << v[(nn++)*n_vars + c] << " ";	    
        }	  
        out << std::endl << std::endl;

#endif

      }
      
    out << "endvars" << std::endl;
  }

  
  // end of the file
  out << std::endl << "endgmv" << std::endl;
}





void GMVIO::add_cell_centered_data (const std::string&       cell_centered_data_name,
				    const std::vector<Real>* cell_centered_data_vals)
{
  libmesh_assert (cell_centered_data_vals != NULL);

  // Make sure there are *at least* enough entries for all the active elements.
  // There can also be entries for inactive elements, they will be ignored.
  libmesh_assert (cell_centered_data_vals->size() >=
	  MeshOutput<MeshBase>::mesh().n_active_elem());
  this->_cell_centered_data[cell_centered_data_name] = cell_centered_data_vals;
}






void GMVIO::read (const std::string& name)
{
  // 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);

  _next_elem_id = 0;

  untested();
  
#ifndef HAVE_GMV

  std::cerr << "Cannot read a GMV file without the GMV API." << std::endl;
  libmesh_error();

#else
  // Clear the mesh so we are sure to start from a pristeen state.
  MeshInput<MeshBase>::mesh().clear();
  
  // It is apparently possible for gmv files to contain