legacy_xdr_io.c

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

// $Id: legacy_xdr_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 <iostream>
#include <iomanip>

#include <vector>
#include <string>

// Local includes
#include "legacy_xdr_io.h"
#include "mesh_base.h"
#include "mesh_data.h"
#include "mesh_tools.h" // MeshTools::n_levels(mesh)
#include "parallel.h"   // - which makes write parallel-only
#include "cell_hex27.h" // Needed for MGF-style Hex27 meshes
#include "boundary_info.h"
#include "libmesh_logging.h"
#include "xdr_mgf.h"
#include "xdr_mesh.h"
#include "xdr_mhead.h"
#include "xdr_soln.h"
#include "xdr_shead.h"

#ifdef USE_COMPLEX_NUMBERS
#include "utility.h"
#endif







// ------------------------------------------------------------
// LegacyXdrIO members
LegacyXdrIO::LegacyXdrIO (MeshBase& mesh, const bool binary) :
  MeshInput<MeshBase> (mesh),
  MeshOutput<MeshBase>(mesh),
  _binary (binary)
{
}




LegacyXdrIO::LegacyXdrIO (const MeshBase& mesh, const bool binary) :
  MeshOutput<MeshBase>(mesh),
  _binary (binary)
{
}




LegacyXdrIO::~LegacyXdrIO ()
{
}




bool & LegacyXdrIO::binary ()
{
  return _binary;
}




bool LegacyXdrIO::binary () const
{
  return _binary;
}



void LegacyXdrIO::read (const std::string& name)
{
  // This is a serial-only process for now;
  // the Mesh should be read on processor 0 and
  // broadcast later
  if (libMesh::processor_id() != 0)
    return;
  
  if (this->binary())
    this->read_binary (name);
  else
    this->read_ascii  (name);
}



void LegacyXdrIO::read_mgf (const std::string& name)
{
  if (this->binary())
    this->read_binary (name, LegacyXdrIO::MGF);
  else
    this->read_ascii  (name, LegacyXdrIO::MGF);
}



void LegacyXdrIO::write (const std::string& name)
{
    if (this->binary())
      this->write_binary (name);
    else
      this->write_ascii  (name);
}



void LegacyXdrIO::write_mgf (const std::string& name)
{
  if (this->binary())
    this->write_binary (name, LegacyXdrIO::MGF);
  else
    this->write_ascii  (name, LegacyXdrIO::MGF);
}



void LegacyXdrIO::read_mgf_soln (const std::string& name,
				 std::vector<Number>& soln,
				 std::vector<std::string>& var_names) const
{
  here();
  std::cerr << "WARNING: this method is deprecated and will disappear soon!"
	    << std::endl;
  
#ifdef USE_COMPLEX_NUMBERS
  
  // buffer for writing separately
  std::vector<Real> real_soln;
  std::vector<Real> imag_soln;
  
  Utility::prepare_complex_data (soln, real_soln, imag_soln);

  this->read_soln (Utility::complex_filename(name, 0), 
		   real_soln, 
		   var_names);
  
  this->read_soln (Utility::complex_filename(name, 1), 
		   imag_soln, 
		   var_names);
  
#else
  
  this->read_soln (name, soln, var_names);
      
#endif
}



void LegacyXdrIO::write_mgf_soln (const std::string& name,
				  std::vector<Number>& soln,
				  std::vector<std::string>& var_names) const
{
  here();
  std::cerr << "WARNING: this method is deprecated and will disappear soon!"
	    << std::endl;
  
#ifdef USE_COMPLEX_NUMBERS
  
  // buffer for writing separately
  std::vector<Real> real_soln;
  std::vector<Real> imag_soln;
  
  Utility::prepare_complex_data (soln, real_soln, imag_soln);
  
  this->write_soln (Utility::complex_filename(name, 0), 
		    real_soln, 
		    var_names);
  
  this->write_soln (Utility::complex_filename(name, 1), 
		    imag_soln, 
		    var_names);
  
#else
  
  this->write_soln (name, soln, var_names);
      
#endif
}



void LegacyXdrIO::read_ascii (const std::string& name, const LegacyXdrIO::FileFormat originator)
{
  // get a writeable reference to the underlying mesh
  MeshBase& mesh = MeshInput<MeshBase>::mesh();
  
  // clear any existing mesh data
  mesh.clear();
    
  // read the mesh
  this->read_mesh (name, originator);
}



void LegacyXdrIO::read_binary (const std::string& name, const LegacyXdrIO::FileFormat originator)
{
#ifndef HAVE_XDR

  std::cerr << "WARNING: Compiled without XDR binary support.\n"
	    << "Will try ASCII instead" << std::endl << std::endl;

  this->read_ascii (name);
  
#else
  
  // get a writeable reference to the underlying mesh
  MeshBase& mesh = MeshInput<MeshBase>::mesh();
  
  // clear any existing mesh data
  mesh.clear();

  // read the mesh
  this->read_mesh (name, originator);
  
#endif
}



void LegacyXdrIO::write_ascii (const std::string& name, const LegacyXdrIO::FileFormat originator)
{
  this->write_mesh (name, originator);
}



void LegacyXdrIO::write_binary (const std::string& name, const LegacyXdrIO::FileFormat originator)
{
#ifndef HAVE_XDR

  std::cerr << "WARNING: Compiled without XDR binary support.\n"
	    << "Will try ASCII instead" << std::endl << std::endl;

  this->write_ascii (name);

#else
  
  this->write_mesh (name, originator);  
  
#endif
}



void LegacyXdrIO::read_mesh (const std::string& name,
			     const LegacyXdrIO::FileFormat originator,
			     MeshData* mesh_data)
{
  // This is a serial-only process for now;
  // the Mesh should be read on processor 0 and
  // broadcast later
  libmesh_assert(libMesh::processor_id() == 0);

  // get a writeable reference to the mesh
  MeshBase& mesh = MeshInput<MeshBase>::mesh();

  // clear any data in the mesh
  mesh.clear();
  
  // Create an XdrMESH object.
  XdrMESH m;

  // Create a pointer
  // to an XdrMESH file
  // header.
  XdrMHEAD mh;

  // Open the XDR file for reading.
  // Note 1: Provide an additional argument
  // to specify the dimension.
  //
  // Note 2: Has to do the right thing for
  // both binary and ASCII files.
  m.set_orig_flag(originator);
  m.init((this->binary() ? XdrMGF::DECODE : XdrMGF::R_ASCII), name.c_str(), 0); // mesh files are always number 0 ...

  // From here on, things depend
  // on whether we are reading or
  // writing!  First, we define
  // header variables that may
  // be read OR written.
  unsigned int              n_blocks = 0;
  unsigned int              n_levels = 0;
  
  if (m.get_orig_flag() == LegacyXdrIO::LIBM)
    n_levels = m.get_num_levels();
  
  
  std::vector<ElemType>     etypes;
  std::vector<unsigned int> neeb;
	
  // Get the information from
  // the header, and place it
  // in the header pointer.
  m.header(&mh);
	
  // Read information from the
  // file header.  This depends on
  // whether its a libMesh or MGF mesh.
  const int numElem     = mh.getNumEl();
  const int numNodes    = mh.getNumNodes();
  const int totalWeight = mh.getSumWghts();
  const int numBCs      = mh.getNumBCs();
  
  // If a libMesh-type mesh, read the augmented mesh information
  if ((m.get_orig_flag() == LegacyXdrIO::DEAL) || (m.get_orig_flag() == LegacyXdrIO::LIBM))
    {
      // Read augmented header
      n_blocks = mh.get_n_blocks();
      
      etypes.resize(n_blocks);
      mh.get_block_elt_types(etypes);
      
      mh.get_num_elem_each_block(neeb); 
    }

  
  
  // Read the connectivity  
  std::vector<int> conn;
  
  // Now that we know the
  // number of nodes and elements,
  // we can resize the
  // appropriate vectors if we are
  // reading information in.
  mesh.reserve_nodes (numNodes);
  mesh.reserve_elem  (numElem);
  
  // Each element stores two extra
  // locations: one which tells
  // what type of element it is,
  // and one which tells how many
  // nodes it has. Therefore,
  // the total number of nodes
  // (totalWeight) must be augmented
  // by 2 times the number of elements
  // in order to read in the entire
  // connectivity array.
  
  // Note: This section now depends on
  // whether we are reading an old-style libMesh, 
  // MGF, or a new-style libMesh mesh.  
  if (m.get_orig_flag() == LegacyXdrIO::DEAL) 
    {
      conn.resize(totalWeight);
      m.Icon(&conn[0], 1, totalWeight);
    }
  
  else if (m.get_orig_flag() == LegacyXdrIO::MGF) 
    {
      conn.resize(totalWeight+(2*numElem));
      m.Icon(&conn[0], 1, totalWeight+(2*numElem));
    }

  else if (m.get_orig_flag() == LegacyXdrIO::LIBM)
    {
      conn.resize(totalWeight);
      m.Icon(&conn[0], 1, totalWeight);
    }
  
  else
    {
      // I don't know what type of mesh it is.
      libmesh_error();
    }

  // read in the nodal coordinates and form points.
  {
    std::vector<Real> coords(numNodes*mesh.spatial_dimension()); // Always use three coords per node
    m.coord(&coords[0], mesh.spatial_dimension(), numNodes);


  
    // Form Nodes out of
    // the coordinates.  If the    
    // MeshData object is active,
    // add the nodes and ids also          
    // to its map.
    for (int innd=0; innd<numNodes; ++innd)      
      {
	Node* node = mesh.add_point (Point(coords[0+innd*3],  
					   coords[1+innd*3],
					   coords[2+innd*3]), innd);
				       
	if (mesh_data != NULL)
	  if (mesh_data->active())
	    {
	      // add the id to the MeshData, so that
	      // it knows the foreign id, even when 
	      // the underlying mesh got re-numbered,
	      // refined, elements/nodes added...   
	      mesh_data->add_foreign_node_id(node, innd);
	    }
      }  
  }

  
  
  // Build the elements.
  // Note: If the originator was MGF, we don't
  // have to do much checking ...
  // all the elements are Hex27.
  // If the originator was
  // this code, we have to loop over
  // et and neeb to read in all the
  // elements correctly.
  //
  // (This used to be before the coords block, but it
  // had to change now that elements store pointers to
  // nodes.  The nodes must exist before we assign them to
  // the elements. BSK, 1/13/2003)
  if ((m.get_orig_flag() == LegacyXdrIO::DEAL) || (m.get_orig_flag() == LegacyXdrIO::LIBM)) 
    {
      unsigned int lastConnIndex = 0;
      unsigned int lastFaceIndex = 0;

      // This map keeps track of elements we've previously
      // constructed, to avoid O(n) lookup times for parent pointers
      // and to enable elements to be added in ascending ID order
      std::map<unsigned int, Elem*> parents;

      // Keep track of Element ids in MGF-style meshes;
      unsigned int _next_elem_id = 0;

      for (unsigned int level=0; level<=n_levels; level++)
      {
        for (unsigned int idx=0; idx<n_blocks; idx++)  
        {
          for (unsigned int e=lastFaceIndex; e<lastFaceIndex+neeb[level*n_blocks+idx]; e++)  
          {
            // Build a temporary element of the right type, so we know how
            // connectivity entries will be on the line for this element.
            AutoPtr<Elem> temp_elem = Elem::build(etypes[idx]);

            // A pointer to the element which will eventually be added to the mesh.
            Elem* elem;

            // New-style libMesh mesh
            if (m.get_orig_flag() == LegacyXdrIO::LIBM)
            {
              unsigned int self_ID   = conn[lastConnIndex + temp_elem->n_nodes()];

#ifdef ENABLE_AMR
              unsigned int parent_ID = conn[lastConnIndex + temp_elem->n_nodes()+1];

              if (level > 0)
              {
                // Do a linear search for the parent
                Elem* my_parent;

                // Search for parent in the parents map (log(n))
                START_LOG("log(n) search for parent", "LegacyXdrIO::read_mesh");
                std::map<unsigned int, Elem*>::iterator it = parents.find(parent_ID);
                STOP_LOG("log(n) search for parent", "LegacyXdrIO::read_mesh");
                
                // If the parent was not previously added, we cannot continue.
                if (it == parents.end())
                {
                  std::cerr << "Parent element with ID " << parent_ID 
                            << " not found." << std::endl; 
                  libmesh_error();
                }

                // Set the my_parent pointer
                my_parent = (*it).second;

                // my_parent is now INACTIVE, since he has children
                my_parent->set_refinement_flag(Elem::INACTIVE);
               
                // Now that we know the parent, build the child
                elem = Elem::build(etypes[idx],my_parent).release();

                // The new child is marked as JUST_REFINED
                elem->set_refinement_flag(Elem::JUST_REFINED); 
                
                // Tell the parent about his new child
                my_parent->add_child(elem);

                // sanity check
                libmesh_assert (my_parent->type() == elem->type());
              }

              // Add level-0 elements to the mesh 
              else
#endif // #ifdef ENABLE_AMR
              {
                elem = Elem::build(etypes[idx]).release();
              }

              // Assign the newly-added element's ID so that future 
              // children which may be added can find it correctly.
              elem->set_id() = self_ID;
                
              // Add this element to the map, it may be a parent for a future element
              START_LOG("insert elem into map", "LegacyXdrIO::read_mesh");
              parents[self_ID] = elem;
              STOP_LOG("insert elem into map", "LegacyXdrIO::read_mesh");
            }

            // MGF-style meshes
            else