unv_io.c

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

// $Id: unv_io.C 2968 2008-08-10 20:38:50Z benkirk $

// 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 <iomanip>
#include <cstdio>   // for std::sprintf
#include <algorithm> // for std::sort
#include <fstream>

// Local includes
#include "libmesh_config.h"
#include "libmesh_logging.h"
#include "unv_io.h"
#include "mesh_data.h"
#include "mesh_base.h"
#include "face_quad4.h"
#include "face_tri3.h"
#include "face_tri6.h"
#include "face_quad8.h"
#include "face_quad9.h"
#include "cell_tet4.h"
#include "cell_hex8.h"
#include "cell_hex20.h"
#include "cell_tet10.h"
#include "cell_prism6.h"

#ifdef HAVE_GZSTREAM
# include "gzstream.h" // For reading/writing compressed streams
#endif





//-----------------------------------------------------------------------------
// UNVIO class static members
const std::string UNVIO::_label_dataset_nodes    = "2411";
const std::string UNVIO::_label_dataset_elements = "2412";



// ------------------------------------------------------------
// UNVIO class members
void UNVIO::clear ()
{
  /*
   * Initialize these to dummy values
   */
  this->_n_nodes     = 0;
  this->_n_elements  = 0;
  this->_need_D_to_e = true;

  this->_assign_nodes.clear();
  this->_ds_position.clear();
}





void UNVIO::read (const std::string& file_name)
{
  if (file_name.rfind(".gz") < file_name.size())
    {
#ifdef HAVE_GZSTREAM
      
      igzstream in_stream (file_name.c_str());
      this->read_implementation (in_stream);
      
#else
      
      std::cerr << "ERROR:  You must have the zlib.h header "
		<< "files and libraries to read and write "
		<< "compressed streams."
		<< std::endl;
      libmesh_error();
      
#endif
      return;      
    }
  
  else
    {
      std::ifstream in_stream (file_name.c_str());
      this->read_implementation (in_stream);
      return;
    }
}


void UNVIO::read_implementation (std::istream& in_stream)
{
  // clear everything, so that
  // we can start from scratch
  this->clear ();


   // Note that we read this file
   // @e twice.  First time to
   // detect the number of nodes
   // and elements (and possible
   // conversion tasks like D_to_e)
   // and the order of datasets
   // (nodes first, then elements,
   // or the other way around),
   // and second to do the actual
   // read.
  std::vector<std::string> order_of_datasets;
  order_of_datasets.reserve(2);
  
  {
    // the first time we read the file,
    // merely to obtain overall info
    if ( !in_stream.good() )
      {
        std::cerr << "ERROR: Input file not good." 
		  << std::endl;
	libmesh_error();
      }

    
    // Count nodes and elements, then let 
    // other methods read the element and 
    // node data.  Also remember which
    // dataset comes first: nodes or elements
    if (this->verbose())
      std::cout << "  Counting nodes and elements" << std::endl;
    

//    bool reached_eof = false;
    bool found_node  = false;
    bool found_elem  = false;


    std::string olds, news;

    while (in_stream.good())
      {
	in_stream >> olds >> news;
	
	// a "-1" followed by a number means the beginning of a dataset
	// stop combing at the end of the file
	while ( ((olds != "-1") || (news == "-1") ) && !in_stream.eof() )
	  {	  
	    olds = news;
	    in_stream >> news;
	  }

//  	if (in_stream.eof())
//  	  {
//  	    reached_eof = true;
//  	    break;
//  	  }

	
	// if beginning of dataset, buffer it in
	// temp_buffer, if desired
	if (news == _label_dataset_nodes)
	  {
	    found_node = true;
	    order_of_datasets.push_back (_label_dataset_nodes);
	    this->count_nodes (in_stream);

	    // we can save some time scanning the file
	    // when we know we already have everything
	    // we want
	    if (found_elem)
	      break;
	  }

	else if (news == _label_dataset_elements)
	  {
	    found_elem = true;
	    order_of_datasets.push_back (_label_dataset_elements);
	    this->count_elements (in_stream);

	    // we can save some time scanning the file
	    // when we know we already have everything
	    // we want
	    if (found_node)
	      break;
	  }
      }


    // Here we should better have found
    // the datasets for nodes and elements,
    // otherwise the unv files is bad!
    if (!found_elem)
      {
	std::cerr << "ERROR: Could not find elements!" << std::endl;
	libmesh_error();
      }

    if (!found_node)
      {
	std::cerr << "ERROR: Could not find nodes!" << std::endl;
	libmesh_error();
      }


    // Don't close, just seek to the beginning
    in_stream.seekg(0, std::ios::beg);
    
    if (!in_stream.good() )
      {
        std::cerr << "ERROR: Cannot re-read input file." 
		  << std::endl;
	libmesh_error();
      }
  }





  // We finished scanning the file,
  // and our member data 
  // \p this->_n_nodes,
  // \p this->_n_elements,
  // \p this->_need_D_to_e
  // should be properly initialized.
  {
    // Read the datasets in the order that
    // we already know
    libmesh_assert (order_of_datasets.size()==2);

    for (unsigned int ds=0; ds < order_of_datasets.size(); ds++)
      {
	if (order_of_datasets[ds] == _label_dataset_nodes)
	  this->node_in    (in_stream);

	else if (order_of_datasets[ds] == _label_dataset_elements)
	  this->element_in (in_stream);

	else
	  libmesh_error();
      }


    // tell the MeshData object that we are finished 
    // reading data
    this->_mesh_data.close_foreign_id_maps ();

    if (this->verbose())
      std::cout << "  Finished." << std::endl << std::endl;
  }
  
  // save memory
  this->_assign_nodes.clear();
  this->_ds_position.clear();
}





void UNVIO::write (const std::string& file_name)
{
  if (file_name.rfind(".gz") < file_name.size())
    {
#ifdef HAVE_GZSTREAM
      
      ogzstream out_stream(file_name.c_str());
      this->write_implementation (out_stream);
      
#else
      
      std::cerr << "ERROR:  You must have the zlib.h header "
		<< "files and libraries to read and write "
		<< "compressed streams."
		<< std::endl;
      libmesh_error();
      
#endif
      
      return;      
    }
  
  else
    {
      std::ofstream out_stream (file_name.c_str());
      this->write_implementation (out_stream);
      return;
    }
}




void UNVIO::write_implementation (std::ostream& out_file)
{
  if ( !out_file.good() )
    {
      std::cerr << "ERROR: Output file not good." 
		<< std::endl;
      libmesh_error();
    }


  MeshBase& mesh = MeshInput<MeshBase>::mesh();
  
  // already know these data, so initialize
  // them.  Does not hurt.
  this->_n_nodes      = mesh.n_nodes();
  this->_n_elements   = mesh.n_elem();
  this->_need_D_to_e  = false;


 
  // we need the MeshData, otherwise we do not
  // know the foreign node id
  if (!this->_mesh_data.active())
    if (!this->_mesh_data.compatibility_mode())
      {
	std::cerr << std::endl
		  << "*************************************************************************" << std::endl
		  << "* WARNING: MeshData neither active nor in compatibility mode.           *" << std::endl
		  << "*          Enable compatibility mode for MeshData.  Use this Universal  *" << std::endl
		  << "*          file with caution: libMesh node and element ids are used.    *" << std::endl
		  << "*************************************************************************" << std::endl
		  << std::endl;
	this->_mesh_data.enable_compatibility_mode();
      }



  // write the nodes,  then the elements
  this->node_out    (out_file);
  this->element_out (out_file);
}





void UNVIO::count_nodes (std::istream& in_file)
{
  START_LOG("count_nodes()","UNVIO");

  // if this->_n_nodes is not 0 the dataset 
  // has already been scanned
  if (this->_n_nodes != 0)
    {
      std::cerr << "Error: Trying to scan nodes twice!" 
		<< std::endl;
      libmesh_error();
    }


  // Read from file, count nodes,
  // check if floats have to be converted
  std::string data;

  in_file >> data; // read the first node label


  if (data == "-1")
    {
      std::cerr << "ERROR: Bad, already reached end of dataset before even starting to read nodes!"
		<< std::endl;
      libmesh_error();
    }

 
  // ignore the misc data for this node
  in_file.ignore(256,'\n');	      


  
  // Now we are there to verify whether we need
  // to convert from D to e or not
  in_file >> data;

  // When this "data" contains a "D", then
  // we have to convert each and every float...
  // But also assume when _this_ specific
  // line does not contain a "D", then the
  // other lines won't, too.
  {
// #ifdef __HP_aCC
//     // Use an "int" instead of unsigned int,
//     // otherwise HP aCC may crash!
//     const int position          = data.find("D",6);
// #else
//     const unsigned int position = data.find("D",6);
// #endif
    std::string::size_type position = data.find("D",6);

    if (position!=std::string::npos) // npos means no position
      {
	this->_need_D_to_e = true;
	
	if (this->verbose())
	  std::cout << "  Convert from \"D\" to \"e\"" << std::endl;
      }
    else
      this->_need_D_to_e = false;
  }

  // read the remaining two coordinates
  in_file >> data;
  in_file >> data;