mesh_data_unv_support.c

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

  if (!_elem_data.empty())
      std::cerr << "WARNING: MeshData currently only supports nodal data for Universal files."
		<< std::endl
		<< "         Will proceed writing only nodal data, ignoring element data."
		<< std::endl;


  /*
   * Write the beginning of the dataset.
   */
  out_file << "    -1\n"
	   << "  " 
	   << _label_dataset_mesh_data
	   << "\n";

  /*
   * Write the header
   */
  if (_unv_header==NULL)
    {
      /*
       * create a header that holds at
       * least sufficient data to specify
       * what this data set currently holds.
       * 
       * The empty constructor automatically
       * takes care of \p dataset_location
       * and \p data_type.
       */
      MeshDataUnvHeader my_header;
      
      /*
       * It remains to set the correct nvaldc...
       */
      my_header.nvaldc = this->n_val_per_node();
      
      /*
       * and the correct data type.  By default
       * only distinguish complex or real data.
       */
#ifdef USE_COMPLEX_NUMBERS
      my_header.data_type = 5;
#else
      my_header.data_type = 2;
#endif

      /*
       * write this default header, then let
       * the AutoPtr go out of scope.  This 
       * will take care of memory management.
       */
      my_header.write (out_file);
    }

  else
    {
      /*
       * make sure our nvaldc coincide.
       */
      if (this->n_val_per_node() != _unv_header->nvaldc)
        {
	  std::cerr << "WARNING: nvaldc=" << _unv_header->nvaldc 
		    << " of attached MeshDataUnvHeader object not valid!" << std::endl
		    << "         re-set nvaldc to " << this->n_val_per_node() << std::endl;
	  _unv_header->nvaldc = this->n_val_per_node();
	}


      /*
       * only issue a warning when data_type does
       * not coincide.  Perhaps user provided some
       * other header in order to convert complex
       * to real...
       */
#ifdef USE_COMPLEX_NUMBERS
      const unsigned int my_data_type = 5;
#else
      const unsigned int my_data_type = 2;
#endif
      if (my_data_type != _unv_header->data_type)
        {
	  std::cerr << "WARNING: data_type=" << _unv_header->data_type 
		    << " of attached MeshDataUnvHeader differs from" << std::endl
		    << "         default value=" << my_data_type
		    << " Perhaps the user wanted this," << std::endl
		    << "         so I use the value from the MeshDataUnvHeader."
		    << std::endl;
	}
      _unv_header->write (out_file);
    }

  
  /*
   * Write the foreign node number and the respective data.
   */
  std::map<const Node*, 
    std::vector<Number> >::const_iterator nit = _node_data.begin();

  char buf[27];
  for (; nit != _node_data.end(); ++nit)
    {
      const Node* node = (*nit).first;

      unsigned int f_n_id = node_to_foreign_id (node);
      std::sprintf(buf, "%10i\n", f_n_id);
      out_file << buf;

      /* since we are iterating over our own map, this libmesh_assert
       * should never break...
       */
      libmesh_assert (this->has_data(node));

      // const reference to the nodal values
      const std::vector<Number>& values = this->get_data(node);

      for (unsigned int v_cnt=0; v_cnt<values.size(); v_cnt++)
	{
#ifdef USE_COMPLEX_NUMBERS
	  std::sprintf(buf, "%13.5E%13.5E", values[v_cnt].real(),
		                            values[v_cnt].imag());
	  out_file << buf;
#else
	  std::sprintf(buf, "%13.5E",
	               static_cast<double>(values[v_cnt]));
	  out_file << buf;
#endif
	}

      out_file << "\n";


    }

  /*
   * Write end of the dataset.
   */
  out_file << "    -1\n";
}





//------------------------------------------------------
// MeshDataUnvHeader functions
MeshDataUnvHeader::MeshDataUnvHeader() :
  dataset_label          (0),
  dataset_name           ("libMesh mesh data"),
  dataset_location       (1),  // default to nodal data
  model_type             (0),          
  analysis_type          (0),
  data_characteristic    (0),
  result_type            (0),
#ifdef USE_COMPLEX_NUMBERS
  data_type              (5),  // default to single precision complex
#else
  data_type              (2),  // default to single precision real
#endif
  nvaldc                 (3),  // default to 3 (principle directions)
  _desired_dataset_label (libMesh::invalid_uint)
{
  id_lines_1_to_5.resize(5);
  std::fill (id_lines_1_to_5.begin(), id_lines_1_to_5.end(), std::string("libMesh default"));
  /*
   * resize analysis specific data.
   */
  record_10.resize(8);
  record_11.resize(2);
  record_12.resize(6);
  record_13.resize(6);
}





MeshDataUnvHeader::~MeshDataUnvHeader()
{
  // empty
}




bool MeshDataUnvHeader::read (std::istream& in_file)
{
  in_file >> this->dataset_label;

  /*
   * currently, we compare only the 
   * dataset_label with the _desired_dataset_label,
   * but it may be easy to also compare the 
   * dataset_name.
   *
   * When the user provided a dataset label, and
   * the current label does _not_ match, then just
   * return false.
   *
   * Otherwise: when the current label matches,
   * or when there is no desired dataset label,
   * simply proceed.
   */
  if ((this->_desired_dataset_label != libMesh::invalid_uint) &&
      (this->dataset_label != this->_desired_dataset_label))
      return false;


  in_file.ignore(256,'\n');
  std::getline(in_file, dataset_name, '\n');

  in_file >> this->dataset_location;
  in_file.ignore(256,'\n');


  for (unsigned int n=0; n<5; n++)
      std::getline(in_file, this->id_lines_1_to_5[n], '\n');


  in_file >> this->model_type     
	  >> this->analysis_type
	  >> this->data_characteristic
	  >> this->result_type
	  >> this->data_type
	  >> this->nvaldc;
    
  for (unsigned int i=0; i<8; i++)
    in_file >> this->record_10[i];

  for (unsigned int i=0; i<2; i++)
    in_file >> this->record_11[i];


  /*
   * There are UNV-files where floats are 
   * written with 'D' as the 10th-power 
   * character. Replace this 'D' by 'e',
   * so that std::atof() can work fine.
   */
  std::string buf;
  in_file >> buf;

  if (need_D_to_e(buf))
    {
      // have to convert _all_ 'D' to 'e'
      this->record_12[0] = std::atof(buf.c_str());

      for (unsigned int i=1; i<6; i++)
        {
	  in_file >> buf;
	  need_D_to_e(buf);
	  this->record_12[i] = std::atof(buf.c_str());
        }

      for (unsigned int i=0; i<6; i++)
        {
	  in_file >> buf;
	  need_D_to_e(buf);
	  this->record_13[i] = std::atof(buf.c_str());
        }
    }
  else
    {
      // no 'D', the stream will recognize the floats
      this->record_12[0] = std::atof(buf.c_str());

      for (unsigned int i=1; i<6; i++)
	  in_file >> this->record_12[i];

      for (unsigned int i=0; i<6; i++)
	  in_file >> this->record_13[i];
    }

  /*
   * no matter whether the user provided a desired
   * dataset label or not: return true, b/c the
   * non-match was already caught before.
   */
  return true;
}




void MeshDataUnvHeader::write (std::ostream& out_file)
{
  
  
  char buf[82];

  std::sprintf(buf, "%6i\n",this->dataset_label);
  
  out_file << buf;
 
  out_file << this->dataset_name << "\n";

  std::sprintf(buf, "%6i\n",this->dataset_location);
  
  out_file << buf;

  for (unsigned int n=0; n<5; n++)
    out_file << this->id_lines_1_to_5[n] << "\n";

  std::sprintf(buf, "%10i%10i%10i%10i%10i%10i\n",
	       model_type,  analysis_type, data_characteristic,
	       result_type, data_type,     nvaldc);
  
  out_file << buf;

  std::sprintf(buf, "%10i%10i%10i%10i%10i%10i%10i%10i\n",
	       record_10[0], record_10[1], record_10[2], record_10[3],
	       record_10[4], record_10[5], record_10[6], record_10[7]);
  
  out_file << buf;

  std::sprintf(buf, "%10i%10i\n", record_11[0], record_11[1]);
  out_file << buf;

  std::sprintf(buf, "%13.5E%13.5E%13.5E%13.5E%13.5E%13.5E\n",
	       static_cast<double>(record_12[0]),
               static_cast<double>(record_12[1]),
               static_cast<double>(record_12[2]),
               static_cast<double>(record_12[3]),
               static_cast<double>(record_12[4]),
               static_cast<double>(record_12[5]));
  
  out_file << buf;

  std::sprintf(buf, "%13.5E%13.5E%13.5E%13.5E%13.5E%13.5E\n",
               static_cast<double>(record_13[0]),
               static_cast<double>(record_13[1]),
               static_cast<double>(record_13[2]),
               static_cast<double>(record_13[3]),
               static_cast<double>(record_13[4]),
               static_cast<double>(record_13[5]));
  
  out_file << buf;
}





bool MeshDataUnvHeader::need_D_to_e (std::string& number)
{
  // find "D" in string, start looking at 6th element, to improve speed.
  // We dont expect a "D" earlier

// #ifdef __HP_aCC
//   // Use an "int" instead of unsigned int,
//   // otherwise HP aCC may crash!
//   const int position = number.find("D",6);
// #else
//   const unsigned int position = number.find("D",6);
// #endif
  std::string::size_type position = number.find("D",6);
    
  if(position!=std::string::npos)     // npos means no position
    {
      // replace "D" in string
      number.replace(position,1,"e"); 
      return true;
    }
  else
      // we assume that if this one number is written correctly, all numbers are
      return false;
}



void MeshDataUnvHeader::which_dataset (const unsigned int ds_label)
{
  this->_desired_dataset_label = ds_label;
}



void MeshDataUnvHeader::operator = (const MeshDataUnvHeader& omduh)
{
  this->dataset_label          = omduh.dataset_label;
  this->dataset_name           = omduh.dataset_name;
  this->dataset_location       = omduh.dataset_location;
  this->id_lines_1_to_5        = omduh.id_lines_1_to_5;

  this->model_type             = omduh.model_type;
  this->analysis_type          = omduh.analysis_type;
  this->data_characteristic    = omduh.data_characteristic;
  this->result_type            = omduh.result_type;

#ifdef USE_COMPLEX_NUMBERS
  /*
   * in complex mode allow only
   * values 5 or 6 (complex) for data_type
   */
  if ((omduh.data_type == 5) ||
      (omduh.data_type == 6))
      this->data_type          = omduh.data_type;
  else
    {
#  ifdef DEBUG
      std::cerr << "WARNING: MeshDataUnvHeader::operator=(): Other object has data_type for" << std::endl
		<< "         real values.  Will use default data_type=5 during assignment." << std::endl
		<< std::endl;
#  endif
      this->data_type          = 5;
    }

#else

  /*
   * in real mode allow only
   * values 2 or 4 (real) for data_type
   */
  if ((omduh.data_type == 2) ||
      (omduh.data_type == 4))
      this->data_type          = omduh.data_type;
  else
    {
#  ifdef DEBUG
      std::cerr << "WARNING: Other MeshDataUnvHeader has data_type for complex values." << std::endl
		<< "         Data import will likely _not_ work and result in infinite loop," << std::endl
		<< "         provided the user forgot to re-size nvaldc to 2*nvaldc_old!" << std::endl
		<< std::endl;
#  endif
      this->data_type          = 2;
    }

#endif

  this->nvaldc                 = omduh.nvaldc;

  this->record_10              = omduh.record_10;
  this->record_11              = omduh.record_11;
  this->record_12              = omduh.record_12;
  this->record_13              = omduh.record_13;

  this->_desired_dataset_label = omduh._desired_dataset_label;
}




bool MeshDataUnvHeader::operator == (const MeshDataUnvHeader& omduh) const
{
  return (this->dataset_label          == omduh.dataset_label       &&
	  this->dataset_name           == omduh.dataset_name        &&
	  this->dataset_location       == omduh.dataset_location    &&
	  this->id_lines_1_to_5        == omduh.id_lines_1_to_5     &&

	  this->model_type             == omduh.model_type          &&
	  this->analysis_type          == omduh.analysis_type       &&
	  this->data_characteristic    == omduh.data_characteristic &&
	  this->result_type            == omduh.result_type         &&

	  this->data_type              == omduh.data_type           &&
	  this->nvaldc                 == omduh.nvaldc              &&

	  this->record_10              == omduh.record_10           &&
	  this->record_11              == omduh.record_11           &&
	  this->record_12              == omduh.record_12           &&
	  this->record_13              == omduh.record_13           &&

	  this->_desired_dataset_label == omduh._desired_dataset_label);
}