plt_loader_write.c

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

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

// Copyright (C) 2002-2007  Benjamin S. Kirk
  
// 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>

// Local includes
#include "plt_loader.h"



// extern "C" {

//   void open_   (const int*, const char*, int*);
//   void idata_  (const int*, int*, int*);
//   void fdata_  (const int*, float*, int*);
//   void ddata_  (const int*, double*, int*);
//   void close_  (const int*);
// }



//-----------------------------------------------------------------------------
// PltLoader write members
// void PltLoader::write_plot3d (const std::string& basename,
// 			      const bool reverse,
// 			      const bool gridonly) const
// {
//   const std::string gname = basename + ".g";
//   const std::string qname = basename + ".q";
  
//   // FORTAN file unit numbers
//   const int gunit = 25;
//   const int qunit = 26;

//   // Tell the user what files we are creating
//   if (this->verbose())
//     {
//       std::cout << "Plot3D output will be written to " << gname;

//       if (!gridonly)
// 	std::cout << " and " << qname;

//       std::cout << std::endl;
//     }

  
//   // Open the FORTRAN unformatted file
//   {
//     libmesh_assert (gname.size() == qname.size());
//     int len = gname.size();
    
//     open_ (&gunit, gname.c_str(), &len);

//     if (!gridonly)
//       open_ (&qunit, qname.c_str(), &len);
//   }
  
//   // Write the headers
//   {
//     std::vector<int> ints;
//     ints.reserve (3*this->n_zones());
    
//     for (unsigned int zn=0; zn<this->n_zones(); zn++)
//       {
// 	ints.push_back(this->imax(zn));
// 	ints.push_back(this->jmax(zn));
// 	ints.push_back(this->kmax(zn));
//       }
    
//     int nb  = this->n_zones();
//     int one = 1;
//     int len = ints.size();

//     libmesh_assert (static_cast<unsigned int>(len) == 3*this->n_zones());
    
//     idata_ (&gunit, &nb, &one);
//     idata_ (&gunit, &ints[0], &len);

//     if (!gridonly)
//       {
// 	idata_ (&qunit, &nb, &one);
// 	idata_ (&qunit, &ints[0], &len);
//       }	
//   }


//   // Variables to write to plot3D file
//   std::vector<unsigned int> write_vars;
//   write_vars.reserve (5);

//   std::fill (write_vars.begin(), write_vars.end(), 0);

  
//   //------------------------------------------------------------------------
//   // Ask the user which variables to write
//   if (!gridonly)
//     {
//       std::cout << "Variables:" << std::endl;
      
//       for (unsigned int v=0; v<this->n_vars(); v++)
// 	std::cout << " " << v << ") \"" << this->var_name(v) << "\""
// 		  << std::endl;
//       std::cout << std::endl;

//       int n_write_vars = 0;

//       while (true)
// 	{
// 	  std::cout << "How many variables to write to the Plot3D file? 1<=n<=" << this->n_vars()
// 		    << " "
// 		    << std::endl
// 		    << "(-1 writes them all): ";
	  
// 	  std::cin >> n_write_vars;

// 	  if (n_write_vars == -1)
// 	    break;

// 	  if ((n_write_vars >= 1) &&
// 	      (n_write_vars <= static_cast<int>(this->n_vars())))
// 	    break;
// 	};


//       // The user wants all the variables
//       if ((n_write_vars == -1) ||
// 	  (n_write_vars == static_cast<int>(this->n_vars())))
// 	{
// 	  for (unsigned int wv=0; wv<this->n_vars(); wv++)
// 	    write_vars.push_back (wv);
// 	}

//       // The user wants a subset of the variables
//       else
// 	{
// 	  libmesh_assert (n_write_vars >= 1);
// 	  libmesh_assert (n_write_vars < static_cast<int>(this->n_vars()));
	  
// 	  std::cout << "Select the " << n_write_vars << " variables to write to the Plot3D file: "
// 		    << std::endl;
	  
// 	  for (int wv=0; wv<n_write_vars; wv++)
// 	    {
// 	      int num=0;
	      
// 	      std::cin >> num;
	      
// 	      libmesh_assert (num < static_cast<int>(this->n_vars()));
	      
// 	      write_vars.push_back (num);
// 	    }
// 	}
      
//       std::cout << std::endl;
//     } // if (!gridonly)


  
//   //------------------------------------------------------------------------
//   // Write the coordinates & data for each block
//   for (unsigned int zn=0; zn<this->n_zones(); zn++)
//     {     
//       // Write the coordinates
//       { 
// 	std::vector<float> coords;   // the nodal coordinates
// 	coords.reserve (3*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
	
// 	for (unsigned int v=0; v<3; v++)
// 	  {
// 	    unsigned int l=0;
	    
// 	    for (unsigned int k=0; k<this->kmax(zn); k++)
// 	      for (unsigned int j=0; j<this->jmax(zn); j++)
// 		for (unsigned int i=0; i<this->imax(zn); i++)
// 		  {
// 		    libmesh_assert (l < _data[zn][v].size());
// 		    coords.push_back (_data[zn][v][l++]);
// 		  }
// 	  }
	
// 	// Write to the grid file
// 	{
// 	  int len = coords.size();
// 	  libmesh_assert (static_cast<unsigned int>(len) ==
// 		  3*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
	  
// 	  fdata_ (&gunit, &coords[0], &len);
// 	}
//       }

      
//       //------------------------------------------------------------------------
//       // Write the data
//       if (!gridonly)
// 	{ 
// 	  std::vector<float> data;     // arbitrary data
// 	  std::vector<float> conds(4); // plot3D conditions [FSMACH, ALPHA, RE, TIME]
// 	  data.reserve (write_vars.size()*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
// 	  std::fill    (conds.begin(), conds.end(), 0.);

// 	  if (zn == 0)
// 	    std::cout << "  Writing ";
	  
// 	  for (unsigned int i=0; i<write_vars.size(); i++)
// 	    {
// 	      // Number of the variable to write
// 	      const unsigned int v = write_vars[i];
	      
// 	      libmesh_assert (v < this->n_vars());

// 	      // Tell the user what variable we are writing, but only
// 	      // once per file.
// 	      if (zn == 0)
// 		std::cout << "\"" << this->var_name(v) << "\" ";
	      
// 	      unsigned int l=0;
	      
// 	      for (unsigned int k=0; k<this->kmax(zn); k++)
// 		for (unsigned int j=0; j<this->jmax(zn); j++)
// 		  for (unsigned int i=0; i<this->imax(zn); i++)
// 		    {
// 		      libmesh_assert (l < _data[zn][v].size());
// 		      data.push_back ((v < this->n_vars()) ?
// 				      _data[zn][v][l++] : 0.);
// 		    }
// 	    }

// 	  if (zn == 0)
// 	    std::cout << "to " << qname << std::endl;
	  
// 	  // Write to the solution file
// 	  {
// 	    int len = conds.size();
	    
// 	    fdata_ (&qunit, &conds[0], &len);
// 	  }
	  
// 	  // Write to the solution file
// 	  {
// 	    int len = data.size();
// 	    libmesh_assert (static_cast<unsigned int>(len) ==
// 		    write_vars.size()*this->imax(zn)*this->jmax(zn)*this->kmax(zn));
	    
// 	    fdata_ (&qunit, &data[0], &len);
// 	  }
// 	}
//     }

//   // Close the FORTAN files
//   close_ (&gunit);

//   if (!gridonly)
//     close_ (&qunit);

//   // Possibly reverse the orders
//   if (reverse)
//     {
//       if (this->verbose())
// 	std::cout << "Reversing byte-ordering for output files."
// 		  << std::endl;
      
//       Utility::reverse_endian (gname);

//       if (!gridonly)
// 	Utility::reverse_endian (qname);
//     }
// }



// void PltLoader::write_tri (const std::string& name,
// 			   const bool reverse,
// 			   const bool gridonly) const
//   {
//   // Check out
//   // http://people.nas.nasa.gov/~aftosmis/cart3d/cart3dTriangulations.html
//   // for the .tri, .triq format
  
//   // FORTRAN file unit numbers
//   const int gunit = 25;

//   if (this->verbose())
//     {
//       std::cout << "Writing unformatted .tri file " << name
// 		<< std::endl;

//       if (gridonly)
// 	std::cout << "Only writing the grid to " << name
// 		  << std::endl;
//     }


//   // Open the FORTRAN unformatted file
//   {
//     int len = name.size();

//     open_ (&gunit, name.c_str(), &len);
//   }


//   // Write the header
//   unsigned int n_nodes =0;
//   unsigned int n_tri   =0;
//   unsigned int n_scalar=this->n_vars()-3;

//   {
//     std::vector<int> ints;
    
//     for (unsigned int zone=0; zone<this->n_zones(); zone++)
//       {
// 	libmesh_assert (this->elem_type(zone) == TRI);
// 	n_nodes += this->n_nodes(zone);
// 	n_tri   += this->n_elem(zone);
//       }
  
//     ints.push_back (n_nodes);
//     ints.push_back (n_tri);

//     if (!gridonly)
//       if (this->n_vars() > 3)
// 	ints.push_back(n_scalar);
  
//     int len = ints.size();
//     idata_ (&gunit, &ints[0], &len);
//   }
  
//   // Write the nodal values.
//   {
//     std::vector<float> coords;
//     coords.reserve (3*n_nodes);

//     for (unsigned int zone=0; zone<this->n_zones(); zone++)
//       for (unsigned int n=0; n<this->n_nodes(zone); n++)
// 	{
// 	  coords.push_back (_data[zone][0][n]);
// 	  coords.push_back (_data[zone][1][n]);
// 	  coords.push_back (_data[zone][2][n]);
// 	}
//     // Got all the nodes for all the zones


//     int len = coords.size();
//     fdata_ (&gunit, &coords[0], &len);
//   }

//   // Write the connectivity
//   {
//     std::vector<int> conn;
//     conn.reserve (3*n_tri);

//     for (unsigned int zone=0; zone<this->n_zones(); zone++)
//       {
// 	// The connectivity for this zone
// 	const std::vector<int> & zconn = _conn[zone];

// 	libmesh_assert (!zconn.empty());
	
// 	// Append the connectivity for this zone to the connectivity
// 	// array
// 	conn.insert (conn.end(), zconn.begin(), zconn.end());
//       }

//     int len = conn.size();
//     libmesh_assert (static_cast<unsigned int>(len) == 3*n_tri);
//     idata_ (&gunit, &conn[0], &len);    
//   }


//   // Write the component index for each triangle
//   {
//     std::vector<int> comp;
//     comp.reserve (n_tri);

//     for (unsigned int zone=0; zone<this->n_zones(); zone++)
//       comp.insert (comp.end(), this->n_elem(zone), zone+1);

//     int len = comp.size();
//     libmesh_assert (static_cast<unsigned int>(len) == n_tri);
//     idata_ (&gunit, &comp[0], &len);
//   }

  
//   // Possibly write additional values for each node
//   if (!gridonly)
//     if (this->n_vars() > 3)
//       {
// 	if (this->verbose())
// 	  {
// 	    std::cout << "Writing variables ";

// 	    for (unsigned int v=3; v<this->n_vars(); v++)
// 	      std::cout << "\"" << this->var_name(v) << "\" ";

// 	    std::cout << "to the output file " << name
// 		      << std::endl;
// 	  }
	
// 	std::vector<float> data;

// 	data.reserve (n_nodes*(this->n_vars()-3));

// 	for (unsigned int zone=0; zone<this->n_zones(); zone++)
// 	  for (unsigned int n=0; n<this->n_nodes(zone); n++)
// 	    for (unsigned int v=3; v<this->n_vars(); v++)
// 	      data.push_back (_data[zone][v][n]);

// 	int len = data.size();
// 	libmesh_assert (static_cast<unsigned int>(len) ==
// 		n_nodes*(this->n_vars()-3));
// 	fdata_ (&gunit, &data[0], &len);
//       }

  
//   // Close the FORTRAN file
//   close_ (&gunit);

  
//   // Possibly reverse the orders
//   if (reverse)
//     {
//       if (this->verbose())
// 	std::cout << "Reversing byte-ordering for output files."
// 		  << std::endl;
      
//       Utility::reverse_endian (name);
//     }
// }



void PltLoader::write_dat (const std::string& name,
			   const unsigned int version) const
{
  std::ofstream out (name.c_str());

  out << "TITLE=\""
      << this->title()
      << "\""
      << '\n';

  out << "VARIABLES = ";

  for (unsigned int v=0; v<this->n_vars(); v++)
    out << "\"" << this->var_name(v) << "\"\n";

  for (unsigned int z=0; z<this->n_zones(); z++)
    {
      out << "ZONE T=\"" << this->zone_name(z) << "\"\n";
      out << " I="  << this->imax(z)
	  << ", J=" << this->jmax(z)
	  << ", K=" << this->kmax(z);

      // Write BLOCK data for this zone
      if (this->zone_type(z) == BLOCK)
	{
	  if (version < 10)
	    {
	      out << ", F=BLOCK\n";
	    }
	  else
	    {
	      out << ", ZONETYPE=Ordered\n"
		  << "DATAPACKING=BLOCK\n";
	    }
      
	  out << "DT=(";
	  for (unsigned int v=0; v<this->n_vars(); v++)
	    out << "SINGLE ";
	  out << ")\n";

	  out.precision(9);
	  
	  for (unsigned int v=0; v<this->n_vars(); v++)
	    {
	      unsigned int l=0;
	      
	      for (unsigned int k=0; k<this->kmax(z); k++)
		for (unsigned int j=0; j<this->jmax(z); j++)
		  for (unsigned int i=0; i<this->imax(z); i++)
		    {
		      // GCC 2.95.3 has scientific in the ios class instead
		      // of in namespace std::
#ifndef BROKEN_IOSTREAM
		      out << std::scientific
			  << _data[z][v][l++] << " ";
#else
		      out << std::ios::scientific
			  << _data[z][v][l++] << " ";
#endif
		      // Throw in a newline every 5 entries to
		      // avoid really long lines.
 		      if (l%5 == 0)
 			out << '\n';
		    }

	      if (l%5 != 0)
		out << '\n';
	    }
	}

      // Write POINT data for this zone
      else if (this->zone_type(z) == POINT)
	{
	  if (version < 10)
	    {
	      out << ", F=POINT\n";
	    }
	  else
	    {
	      out << ", ZONETYPE=Ordered\n"
		  << "DATAPACKING=POINT\n";
	    }
	  
	  out << "DT=(";
	  for (unsigned int v=0; v<this->n_vars(); v++)
	    out << "SINGLE ";
	  out << ")\n";
	  
	  out.precision(9);
	  
	  {
	    unsigned int l=0;
	    
	    for (unsigned int k=0; k<this->kmax(z); k++)
	      for (unsigned int j=0; j<this->jmax(z); j++)
		for (unsigned int i=0; i<this->imax(z); i++)
		  {
		    for (unsigned int v=0; v<this->n_vars(); v++)

		      // GCC 2.95.3 has scientific in the ios class instead
		      // of in namespace std::
#ifndef BROKEN_IOSTREAM
		      out << std::scientific
			  << _data[z][v][l] << " ";
#else
		      out << std::ios::scientific
			  << _data[z][v][l] << " ";
#endif
		    out << '\n';
		    
		    l++;
		  }
	  }
	}
    }
}