gmv_io.c

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

#ifdef USE_COMPLEX_NUMBERS

          // in case of complex data, write _three_ 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";

#endif
}






void GMVIO::write_ascii_old_impl (const std::string& fname,
				  const std::vector<Number>* v,
				  const std::vector<std::string>* solution_names)
{
  // Open the output file stream
  std::ofstream out (fname.c_str());
  
  libmesh_assert (out.good());

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

  // Make sure our nodes are contiguous and serialized
  libmesh_assert (mesh.n_nodes() == mesh.max_node_id());
  // libmesh_assert (mesh.is_serial());
  if (!mesh.is_serial())
    {
      if (libMesh::processor_id() == 0)
        std::cerr << "Error: GMVIO cannot yet write a ParallelMesh solution"
                  << std::endl;
      return;
    }

  unsigned int mesh_max_p_level = 0;
  
  // Begin interfacing with the GMV data file

  // FIXME - if subdivide_second_order() is off,
  // we probably should only be writing the
  // vertex nodes - RHS
  {
    // write the nodes
    
    out << "gmvinput ascii\n\n";
    out << "nodes " << mesh.n_nodes() << '\n';
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.point(v)(0) << " ";
         
    out << '\n';
    
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.point(v)(1) << " ";
    
    out << '\n';
    
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.point(v)(2) << " ";
     
    out << '\n' << '\n';
  }


  
  {
    // write the connectivity
    
    out << "cells ";
    if (this->subdivide_second_order())
      out << mesh.n_active_sub_elem();
    else
      out << mesh.n_active_elem();
    out << '\n';

    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:
	{
	  // The same temporary storage will be used for each element
	  std::vector<unsigned int> conn;

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

              if (this->subdivide_second_order())
	        for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
	          {
		    out << "line 2\n";
		    (*it)->connectivity(se, TECPLOT, conn);
		    for (unsigned int i=0; i<conn.size(); i++)
		      out << conn[i] << " ";
		
		    out << '\n';
	          }
              else
                {
		  out << "line 2\n";
                  if ((*it)->default_order() == FIRST)
		    (*it)->connectivity(0, TECPLOT, conn);
                  else
                    {
                      AutoPtr<Elem> lo_elem = Elem::build(
                        Elem::first_order_equivalent_type((*it)->type()));
                      for (unsigned int i = 0; i != lo_elem->n_nodes(); ++i)
                        lo_elem->set_node(i) = (*it)->get_node(i);
		      lo_elem->connectivity(0, TECPLOT, conn);
                    }
		  for (unsigned int i=0; i<conn.size(); i++)
		    out << conn[i] << " ";
		
		  out << '\n';
                }
            }
	  break;
	}
	
      case 2:
	{
	  // The same temporary storage will be used for each element
	  std::vector<unsigned int> conn;
	  
	  for ( ; it != end; ++it)
            {
              mesh_max_p_level = std::max(mesh_max_p_level,
                                          (*it)->p_level());

              if (this->subdivide_second_order())
	        for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
	          {
		    // Quad elements
		    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\n";
		        (*it)->connectivity(se, TECPLOT, conn);
		        for (unsigned int i=0; i<conn.size(); i++)
		          out << conn[i] << " ";
		      }

		    // Triangle elements
		    else if (((*it)->type() == TRI3) ||
			     ((*it)->type() == TRI6))
		      {
		        out << "tri 3\n";
		        (*it)->connectivity(se, TECPLOT, conn);
		        for (unsigned int i=0; i<3; i++)
		          out << conn[i] << " ";
		      }
		    else
		      libmesh_error();
                  }
              else // !this->subdivide_second_order()
                {
		  // Quad elements
		  if (((*it)->type() == QUAD4)
#ifdef ENABLE_INFINITE_ELEMENTS
		      || ((*it)->type() == INFQUAD4)
#endif
		      )
                    {
		      (*it)->connectivity(0, TECPLOT, conn);
		      out << "quad 4\n";
		      for (unsigned int i=0; i<conn.size(); i++)
		        out << conn[i] << " ";
		    }
		  else if (((*it)->type() == QUAD8) ||
		           ((*it)->type() == QUAD9)
#ifdef ENABLE_INFINITE_ELEMENTS
		           || ((*it)->type() == INFQUAD6)
#endif
		          )
		    {
                      AutoPtr<Elem> lo_elem = Elem::build(
                        Elem::first_order_equivalent_type((*it)->type()));
                      for (unsigned int i = 0; i != lo_elem->n_nodes(); ++i)
                        lo_elem->set_node(i) = (*it)->get_node(i);
		      lo_elem->connectivity(0, TECPLOT, conn);
		      out << "quad 4\n";
		      for (unsigned int i=0; i<conn.size(); i++)
		        out << conn[i] << " ";
		    }
		  else if ((*it)->type() == TRI3)
		    {
		      (*it)->connectivity(0, TECPLOT, conn);
		      out << "tri 3\n";
		      for (unsigned int i=0; i<3; i++)
		        out << conn[i] << " ";
		    }
		  else if ((*it)->type() == TRI6)
		    {
                      AutoPtr<Elem> lo_elem = Elem::build(
                        Elem::first_order_equivalent_type((*it)->type()));
                      for (unsigned int i = 0; i != lo_elem->n_nodes(); ++i)
                        lo_elem->set_node(i) = (*it)->get_node(i);
		      lo_elem->connectivity(0, TECPLOT, conn);
		      out << "tri 3\n";
		      for (unsigned int i=0; i<3; i++)
		        out << conn[i] << " ";
		    }
		
		  out << '\n';
	        }
            }
	  
	  break;
	}
	
	
      case 3:
	{
	  // The same temporary storage will be used for each element
	  std::vector<unsigned int> conn;
	  
	  for ( ; it != end; ++it)
            {
              mesh_max_p_level = std::max(mesh_max_p_level,
                                          (*it)->p_level());

              if (this->subdivide_second_order())
	        for (unsigned int se=0; se<(*it)->n_sub_elem(); se++)
	          {

#ifndef  ENABLE_INFINITE_ELEMENTS
		    if (((*it)->type() == HEX8)   ||    
		        ((*it)->type() == HEX27))
		      {
		        out << "phex8 8\n";
		        (*it)->connectivity(se, TECPLOT, conn);
		        for (unsigned int i=0; i<conn.size(); i++)
		          out << conn[i] << " ";
		      }
		
		    else if ((*it)->type() == HEX20)
		      {
		        out << "phex20 20\n";
		        out << (*it)->node(0)+1  << " "
			    << (*it)->node(1)+1  << " "
			    << (*it)->node(2)+1  << " "
			    << (*it)->node(3)+1  << " "
			    << (*it)->node(4)+1  << " "
			    << (*it)->node(5)+1  << " "
			    << (*it)->node(6)+1  << " "
			    << (*it)->node(7)+1  << " "
			    << (*it)->node(8)+1  << " "
			    << (*it)->node(9)+1  << " "
			    << (*it)->node(10)+1 << " "
			    << (*it)->node(11)+1 << " "
			    << (*it)->node(16)+1 << " "
			    << (*it)->node(17)+1 << " "
			    << (*it)->node(18)+1 << " "
			    << (*it)->node(19)+1 << " "
			    << (*it)->node(12)+1 << " "
			    << (*it)->node(13)+1 << " "
			    << (*it)->node(14)+1 << " "
			    << (*it)->node(15)+1 << " ";
		      }
#else
		    /*
		     * In case of infinite elements, HEX20
		     * should be handled just like the
		     * INFHEX16, since these connect to each other
		     */
		    if (((*it)->type() == HEX8)     ||
		        ((*it)->type() == HEX27)    ||
		        ((*it)->type() == INFHEX8)  ||
		        ((*it)->type() == INFHEX16) ||
		        ((*it)->type() == INFHEX18) ||
		        ((*it)->type() == HEX20))
		      {
		        out << "phex8 8\n";
		        (*it)->connectivity(se, TECPLOT, conn);
		        for (unsigned int i=0; i<conn.size(); i++)
		          out << conn[i] << " ";
		      }
#endif
		
		    else if (((*it)->type() == TET4)  ||
			     ((*it)->type() == TET10))
		      {
		        out << "tet 4\n";
		        (*it)->connectivity(se, TECPLOT, conn);
		        out << conn[0] << " "
			    << conn[2] << " "
			    << conn[1] << " "
			    << conn[4] << " ";
		      }
#ifndef  ENABLE_INFINITE_ELEMENTS
		    else if (((*it)->type() == PRISM6)  ||
			     ((*it)->type() == PRISM15) ||
			     ((*it)->type() == PRISM18))
#else
		    else if (((*it)->type() == PRISM6)     ||
			     ((*it)->type() == PRISM15)    ||
			     ((*it)->type() == PRISM18)    ||
			     ((*it)->type() == INFPRISM6)  ||
			     ((*it)->type() == INFPRISM12))
#endif
		      {
		        /**
		         * Note that the prisms are treated as
		         * degenerated phex8's.
		         */
		        out << "phex8 8\n";
		        (*it)->connectivity(se, TECPLOT, conn);
		        for (unsigned int i=0; i<conn.size(); i++)
		          out << conn[i] << " ";
		      }
		
		    else
		      {
		        std::cout << "Encountered an unrecognized element "
			          << "type: " << (*it)->type()
				  << "\nPossibly a dim-1 dimensional "
			          << "element?  Aborting..."
			          << std::endl;
		        libmesh_error();
		      }
		
		    out << '\n';
	          }
              else // !this->subdivide_second_order()
                {
                  AutoPtr<Elem> lo_elem = Elem::build(
                    Elem::first_order_equivalent_type((*it)->type()));
                  for (unsigned int i = 0; i != lo_elem->n_nodes(); ++i)
                    lo_elem->set_node(i) = (*it)->get_node(i);
		  if ((lo_elem->type() == HEX8)
#ifdef  ENABLE_INFINITE_ELEMENTS
                      || (lo_elem->type() == HEX27)
#endif
                     )
		    {
		      out << "phex8 8\n";
		      lo_elem->connectivity(0, TECPLOT, conn);
		      for (unsigned int i=0; i<conn.size(); i++)
		        out << conn[i] << " ";
		    }
		
		  else if (lo_elem->type() == TET4)
		    {
		      out << "tet 4\n";
		      lo_elem->connectivity(0, TECPLOT, conn);
		      out << conn[0] << " "
			  << conn[2] << " "
			  << conn[1] << " "
			  << conn[4] << " ";
		    }
		  else if ((lo_elem->type() == PRISM6)
#ifdef  ENABLE_INFINITE_ELEMENTS
                           || (lo_elem->type() == INFPRISM6)
#endif
			   )
		    {
		      /**
		       * Note that the prisms are treated as
		       * degenerated phex8's.
		       */
		      out << "phex8 8\n";
		      lo_elem->connectivity(0, TECPLOT, conn);
		      for (unsigned int i=0; i<conn.size(); i++)
		        out << conn[i] << " ";
		    }
		
		  else
		    {
		      std::cout << "Encountered an unrecognized element "
			        << "type.  Possibly a dim-1 dimensional "
			        << "element?  Aborting..."
			        << std::endl;
		      libmesh_error();
		    }
		
		  out << '\n';
	        }
            }
	  
	  break;
	}
	
      default:
	libmesh_error();
      }
    
    out << '\n';
  }
  

  
  // optionally write the partition information
  if (this->partitioning())
    {
      out << "material "
	  << mesh.n_partitions()
	  << " 0"<< '\n';

      for (unsigned int proc=0; proc<mesh.n_partitions(); proc++)
	out << "proc_" << proc << '\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)->processor_id()+1 << '\n';