mesh_data_unv_support.c

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

// $Id: mesh_data_unv_support.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 <cstdio> // for std::sprintf
#include <fstream>

// Local includes
#include "libmesh_config.h"
#include "mesh_data.h"
#include "auto_ptr.h"

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


//------------------------------------------------------
// MeshData UNV support functions
void MeshData::read_unv (const std::string& file_name)
{
  /*
   * we should better be active or in compatibility mode
   */
  libmesh_assert (this->_active || this->_compatibility_mode);

  /*
   * When reading data, make sure the id maps are ok
   */
  libmesh_assert (this->_node_id_map_closed);
  libmesh_assert (this->_elem_id_map_closed);

  /*
   * clear the data, but keep the id maps
   */
  this->clear();

  /*
   * We can read either ".unv", or ".unv.gz"
   * files, provided zlib.h is there
   */
  if (file_name.rfind(".gz") < file_name.size())
    {
#ifdef HAVE_GZSTREAM
      igzstream in_stream(file_name.c_str());
      this->read_unv_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_unv_implementation (in_stream);
      return;
    }
}






void MeshData::read_unv_implementation (std::istream& in_file)
{
  /*
   * This is the actual implementation of
   * reading in UNV format.  This enables
   * to read either through the conventional
   * C++ stream, or through a stream that
   * allows to read .gz'ed files.
   */
  if ( !in_file.good() )
    {
      std::cerr << "ERROR: Input file not good." 
		<< std::endl;
      libmesh_error();
    }

  const std::string _label_dataset_mesh_data = "2414";

  /*
   * locate the beginning of data set
   * and read it.
   */
  {
    std::string olds, news;

    while (true)
    {
      in_file >> olds >> news;

      /*
       * Yes, really dirty:
       *
       * When we found a dataset, and the user does
       * not want this dataset, we jump back here
       */
    go_and_find_the_next_dataset:

      /*
       * 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_file.eof() )
	{
	  olds = news;
	  in_file >> news;
	}

      if(in_file.eof())
	break;

      /*
       * if beginning of dataset
       */
      if (news == _label_dataset_mesh_data)
        {

	  /*
	   * Now read the data of interest.
	   * Start with the header.  For 
	   * explanation of the variable
	   * dataset_location, see below.
	   */
	  unsigned int dataset_location;
	      
	  /*
	   * the type of data (complex, real,
	   * float, double etc, see below)
	   */
	  unsigned int data_type;

	  /*
	   * the number of floating-point values per entity
	   */
	  unsigned int NVALDC;


	  /*
	   * If there is no MeshDataUnvHeader object
	   * attached
	   */
	  if (_unv_header==NULL)
	    {
	      /*
	       * Ignore the first lines that stand for
	       * analysis dataset label and name.
	       */ 
	      for(unsigned int i=0; i<3; i++)
		in_file.ignore(256,'\n');	      
	      
	      /*
	       * Read the dataset location, where
	       * 1: Data at nodes
	       * 2: Data on elements
	       * other sets are currently not supported.
	       */
	      in_file >> dataset_location;

	      /*
	       * Ignore five ID lines.
	       */ 
	      for(unsigned int i=0; i<6; i++)
		in_file.ignore(256,'\n');
	      
	      /*
	       * These data are all of no interest to us...
	       */
	      unsigned int model_type,          
		  analysis_type,
		  data_characteristic,
		  result_type;

	      /*
	       * Read record 9.
	       */
	      in_file >> model_type           // not used here
		      >> analysis_type        // not used here
		      >> data_characteristic  // not used here
		      >> result_type          // not used here
		      >> data_type
		      >> NVALDC;
	      
	      
	      /*
	       * Ignore record 10 and 11
	       * (Integer analysis type specific data).
	       */ 
	      for (unsigned int i=0; i<3; i++)
		in_file.ignore(256,'\n');
	      
	      /*
	       * Ignore record 12 and record 13.  Since there
	       * exist UNV files with 'D' instead of 'e' as 
	       * 10th-power char, it is safer to use a string
	       * to read the dummy reals.
	       */
	      {
	        std::string dummy_Real;
		for (unsigned int i=0; i<12; i++)
		    in_file >> dummy_Real;
	      }

	    }
	  else
	    {

	      /*
	       * the read() method returns false when
	       * the user wanted a special header, and
	       * when the current header is _not_ the correct
	       * header
	       */
	      if (_unv_header->read(in_file))
	        {
		  dataset_location = _unv_header->dataset_location;
		  NVALDC = _unv_header->nvaldc;
		  data_type = _unv_header->data_type;
		}
	      else
	        {
		  /*
		   * This is not the correct header.  Go
		   * and find the next.  For this to
		   * work correctly, shift to the 
		   * next line, so that the "-1"
		   * disappears from olds 
		   */
		  olds = news;
		  in_file >> news;

		  /*
		   * No good style, i know...
		   */
		  goto go_and_find_the_next_dataset;
		}

	    }

	  /*
	   * Check the location of the dataset.
	   */
	  if (dataset_location != 1)
	    {
	      std::cerr << "ERROR: Currently only Data at nodes is supported." 
			<< std::endl;
	      libmesh_error();
	    }


	  /*
	   * Now get the foreign node id number and the respective nodal data.
	   */
	  int f_n_id;
	  std::vector<Number> values;

	  while(true)
	    {
	      in_file >> f_n_id;
	  
	      /*
	       * if node_nr = -1 then we have reached the end of the dataset.
	       */
	      if (f_n_id==-1)
		  break;

	      /*
	       * Resize the values vector (usually data in three
	       * principle directions, i.e. NVALDC = 3).
	       */
	      values.resize(NVALDC);
	  
	      /*
	       * Read the meshdata for the respective node.
	       */	       
	      for (unsigned int data_cnt=0; data_cnt<NVALDC; data_cnt++)
		{
		  /*
		   * Check what data type we are reading.
		   * 2,4: Real
		   * 5,6: Complex
		   * other data types are not supported yet.
		   * As again, these floats may also be written
		   * using a 'D' instead of an 'e'.
		   */
		  if (data_type == 2 || data_type == 4)
		    {
		      std::string buf;
		      in_file >> buf;
		      MeshDataUnvHeader::need_D_to_e(buf);
#ifdef USE_COMPLEX_NUMBERS
		      values[data_cnt] = Complex(std::atof(buf.c_str()), 0.);
#else
		      values[data_cnt] = std::atof(buf.c_str());
#endif
		    }

		  else if(data_type == 5 || data_type == 6)

		    {
#ifdef USE_COMPLEX_NUMBERS
		      Real re_val, im_val;

		      std::string buf;
		      in_file >> buf;

		      if (MeshDataUnvHeader::need_D_to_e(buf))
		        {
			  re_val = std::atof(buf.c_str());
			  in_file >> buf;
			  MeshDataUnvHeader::need_D_to_e(buf);
			  im_val = std::atof(buf.c_str()); 
			}
		      else
		        {
			  re_val = std::atof(buf.c_str());
			  in_file >> im_val;
			}

		      values[data_cnt] = Complex(re_val,im_val);
#else

		      std::cerr << "ERROR: Complex data only supported" << std::endl
				<< "when libMesh is configured with --enable-complex!"
				<< std::endl;
		      libmesh_error();
#endif
		    }

		  else
		    {
		      std::cerr << "ERROR: Data type not supported." 
				<< std::endl;
		      libmesh_error();
		    }

		} // end loop data_cnt

	      /*
	       * Add the values vector to the MeshData data structure.
	       */
	      const Node* node = foreign_id_to_node(f_n_id);
	      _node_data.insert (std::make_pair(node, values));

	    } // while(true)
	}


      else
        {
	  /*
	   * all other datasets are ignored
	   */
        }

    }
  }


  /*
   * finished reading.  Ready for use, provided
   * there was any data contained in the file.
   */
  libmesh_assert ((this->_node_data.size() != 0) || (this->_elem_data.size() != 0));

  this->_node_data_closed = true;
  this->_elem_data_closed = true;
}






void MeshData::write_unv (const std::string& file_name)
{
  /*
   * we should better be active or in compatibility mode
   */
  libmesh_assert (this->_active || this->_compatibility_mode);

  /*
   * make sure the id maps are ready
   * and that we have data to write
   */
  libmesh_assert (this->_node_id_map_closed);
  libmesh_assert (this->_elem_id_map_closed);

  libmesh_assert (this->_node_data_closed);
  libmesh_assert (this->_elem_data_closed);

  if (file_name.rfind(".gz") < file_name.size())
    {
#ifdef HAVE_GZSTREAM
      ogzstream out_stream(file_name.c_str());
      this->write_unv_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_unv_implementation (out_stream);
      return;
    }
}






void MeshData::write_unv_implementation (std::ostream& out_file)
{
  /*
   * This is the actual implementation of writing
   * unv files, either as .unv or as .unv.gz file
   */
  if ( !out_file.good() )
    {
      std::cerr << "ERROR: Output file not good." 
		<< std::endl;
      libmesh_error();
    }


  /*
   * the beginning marker of the dataset block for
   * nodal/element-associated data (not to be confused
   * with _desired_dataset_label!)
   */
  const std::string  _label_dataset_mesh_data    = "2414";

  /*
   * Currently this function handles only nodal data.
   */
  libmesh_assert (!_node_data.empty());