gnuplot_io.c

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

// $Id: gnuplot_io.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

// C++ includes
#include <fstream>
#include <sstream>
#include <map>

// Local includes
#include "elem.h"
#include "mesh_base.h"
#include "gnuplot_io.h"

GnuPlotIO::GnuPlotIO(const MeshBase& mesh,
		     const std::string& title, 
                     int mesh_properties)
  :
  MeshOutput<MeshBase> (mesh),
  _title(title)
{
  _grid       = (mesh_properties & GRID_ON);
  _png_output = (mesh_properties & PNG_OUTPUT);
}

void GnuPlotIO::write(const std::string& fname)
{

  if(libMesh::processor_id() == 0)
    this->write_solution(fname);

}

void GnuPlotIO::write_nodal_data (const std::string& fname,
				  const std::vector<Number>& soln,
				  const std::vector<std::string>& names)
{
  if(libMesh::processor_id() == 0)
    this->write_solution(fname, &soln, &names);

}




void GnuPlotIO::write_solution(const std::string& fname,
   			       const std::vector<Number>* soln,
			       const std::vector<std::string>* names)
{
  libmesh_assert(libMesh::processor_id() == 0);

  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();

  std::stringstream data_stream_name;
  data_stream_name << fname << "_data";
  const std::string data_file_name = data_stream_name.str();
  
  // This class is designed only for use with 1D meshes
  libmesh_assert (mesh.mesh_dimension() == 1);

  // Make sure we have a solution to plot
  libmesh_assert ((names != NULL) && (soln != NULL));

  // Create an output stream for script file
  std::ofstream out(fname.c_str());

  if (!out.good())
    {
      std::cerr << "ERROR: opening output file " << fname
		<< std::endl;
      libmesh_error();
    }

  // The number of variables in the equation system
  const unsigned int n_vars = names->size();
  
  // Write header to stream
  out << "# This file was generated by gnuplot_io.C\n"
      << "# Stores 1D solution data in GNUplot format\n"
      << "# Execute this by loading gnuplot and typing "
      << "\"call '" << fname << "'\"\n"
      << "reset\n"
      << "set title \"" << _title << "\"\n"
      << "set xlabel \"x\"\n"
      << "set autoscale\n"
      << "set xtics nomirror\n";
      
  if(_grid)
  {
    // construct string for xtic positions at element edges
    std::string xtics;
    std::stringstream xtics_stream;

    MeshBase::const_element_iterator it = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_it =
      mesh.active_local_elements_end();

    unsigned int count = 0;
    unsigned int n_active_elem = mesh.n_active_elem();

    for( ; it != end_it; ++it)
    {
      Elem* el = *it;
      
      // if el is the left edge of the mesh, print its left node position
      if(el->neighbor(0) == NULL)
      {
        xtics_stream << "\"\" " << (*(el->get_node(0)))(0) << ", \\\n";
      }
      xtics_stream << "\"\" " << (*(el->get_node(1)))(0);

      if(count+1 != n_active_elem)
      {
         xtics_stream << ", \\\n";
      }
      count++;
    }
    xtics = xtics_stream.str();
    
  
    out << "set x2tics (" << xtics << ")\nset grid noxtics noytics x2tics\n";
  }

  if(_png_output)
  {
    out << "set terminal png\n";
    out << "set output \"" << fname << ".png\"\n";
  }

  out << "plot "
      << axes_limits
      << " \"" << data_file_name << "\" using 1:2 title \"" << (*names)[0]
      << "\" with lines";
  if(n_vars > 1)
  {
    for(unsigned int i=1; i<n_vars; i++)
    {
      out << ", \\\n\"" << data_file_name << "\" using 1:" << i+2 
          << " title \"" << (*names)[i] << "\" with lines";
    }
  }

  out.close();


  // Create an output stream for data file
  std::ofstream data(data_file_name.c_str());

  if (!data.good())
    {
      std::cerr << "ERROR: opening output data file " << std::endl;
      libmesh_error();
    }

  data << "# This file contains libMesh solution data "
       << "to be plotted in GNUplot\n\n";

  // get ordered nodal data using a map
  typedef std::pair<Real, std::vector<Number> > key_value_pair;
  typedef std::map<Real, std::vector<Number> > map_type;
  typedef map_type::iterator map_iterator;

  map_type node_map;


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

  for ( ; it != end; ++it)
  {
    const Elem* elem = *it;

    for(unsigned int i=0; i<elem->n_nodes(); i++)
    {
      std::vector<Number> values;

      // Get the global id of the node
      unsigned int global_id = elem->node(i);

      for(unsigned int c=0; c<n_vars; c++)
      {
        values.push_back( (*soln)[global_id*n_vars + c] );
      }

      node_map[ mesh.point(global_id)(0) ] = values;
    }
  }

    
  map_iterator map_it = node_map.begin();
  const map_iterator end_map_it = node_map.end();

  for( ; map_it != end_map_it; ++map_it)
  {
    key_value_pair kvp = *map_it;
    std::vector<Number> values = kvp.second;

    data << kvp.first << "\t";

    for(unsigned int i=0; i<values.size(); i++)
    {
      data << values[i] << "\t";
    }

    data << "\n";
  }

  data.close();
}