Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

equation_systems_io.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// 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


#include "libmesh/libmesh_common.h"
#include "libmesh/libmesh_logging.h"


// C++ Includes
#include <cstdio> // for std::sprintf
#include <sstream>

// Local Includes
#include "libmesh/libmesh_version.h"
#include "libmesh/equation_systems.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/parallel_mesh.h"
#include "libmesh/parallel.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/xdr_cxx.h"
#include "libmesh/mesh_refinement.h"

namespace libMesh
{

// Forward Declarations

// Anonymous namespace for implementation details.
namespace {
  std::string local_file_name (const std::string &name)
  {
    std::string basename(name);
    char buf[256];

    if (basename.size() - basename.rfind(".bz2") == 4)
      {
        basename.erase(basename.end()-4, basename.end());
        std::sprintf(buf, "%s.%04d.bz2", basename.c_str(), libMesh::processor_id());
      }
    else if (basename.size() - basename.rfind(".gz") == 3)
      {
        basename.erase(basename.end()-3, basename.end());
        std::sprintf(buf, "%s.%04d.gz", basename.c_str(), libMesh::processor_id());
      }
    else
      std::sprintf(buf, "%s.%04d", basename.c_str(), libMesh::processor_id());

    return std::string(buf);
  }
}




// ------------------------------------------------------------
// EquationSystem class implementation
void EquationSystems::read (const std::string& name,
                            const unsigned int read_flags)
{
  libMeshEnums::XdrMODE mode = READ;
  if (name.find(".xdr") != std::string::npos)
    mode = DECODE;
  this->read(name, mode, read_flags);

#ifdef LIBMESH_ENABLE_AMR
  MeshRefinement mesh_refine(_mesh);
  mesh_refine.clean_refinement_flags();
#endif
}



void EquationSystems::read (const std::string& name,
                            const libMeshEnums::XdrMODE mode,
                            const unsigned int read_flags)
{
#ifdef LIBMESH_ENABLE_EXCEPTIONS

  // If we have exceptions enabled we can be considerate and try
  // to read old restart files which contain infinite element
  // information but do not have the " with infinite elements"
  // string in the version information.

  // First try the read the user requested
  try
    {
      this->_read_impl (name, mode, read_flags);
    }

  // If that fails, try it again but explicitly request we look for infinite element info
  catch (...)
    {
      libMesh::out << "\n*********************************************************************\n"
                    << "READING THE FILE \"" << name << "\" FAILED.\n"
                    << "It is possible this file contains infinite element information,\n"
                    << "but the version string does not contain \" with infinite elements\"\n"
                    << "Let's try this again, but looking for infinite element information...\n"
                    << "*********************************************************************\n"
                    << std::endl;

      try
        {
          this->_read_impl (name, mode, read_flags | EquationSystems::TRY_READ_IFEMS);
        }

      // If all that failed, we are out of ideas here...
      catch (...)
        {
          libMesh::out << "\n*********************************************************************\n"
                        << "Well, at least we tried!\n"
                        << "Good Luck!!\n"
                        << "*********************************************************************\n"
                        << std::endl;
          throw;
        }
    }

#else

  // no exceptions - cross your fingers...
  this->_read_impl (name, mode, read_flags);

#endif // #ifdef LIBMESH_ENABLE_EXCEPTIONS

#ifdef LIBMESH_ENABLE_AMR
  MeshRefinement mesh_refine(_mesh);
  mesh_refine.clean_refinement_flags();
#endif
}



void EquationSystems::_read_impl (const std::string& name,
                                  const libMeshEnums::XdrMODE mode,
                                  const unsigned int read_flags)
{
   // Set booleans from the read_flags argument
   const bool read_header          = read_flags & EquationSystems::READ_HEADER;
   const bool read_data            = read_flags & EquationSystems::READ_DATA;
   const bool read_additional_data = read_flags & EquationSystems::READ_ADDITIONAL_DATA;
   const bool read_legacy_format   = read_flags & EquationSystems::READ_LEGACY_FORMAT;
   const bool try_read_ifems       = read_flags & EquationSystems::TRY_READ_IFEMS;
   const bool read_basic_only      = read_flags & EquationSystems::READ_BASIC_ONLY;
         bool read_parallel_files  = false;

  std::map<std::string, System*> xda_systems;

  // This will unzip a file with .bz2 as the extension, otherwise it
  // simply returns the name if the file need not be unzipped.
  Xdr io ((libMesh::processor_id() == 0) ? name : "", mode);
  libmesh_assert (io.reading());

  {
    // 1.)
    // Read the version header.
    std::string version = "legacy";
    if (!read_legacy_format)
      {
        if (libMesh::processor_id() == 0) io.data(version);
        CommWorld.broadcast(version);

        // All processors have the version header, if it does not contain
        // "libMesh" something then it is a legacy file.
        std::string::size_type lm_pos = version.find("libMesh");
        if (!(lm_pos < version.size()))
          {
            io.close();

            // Recursively call this read() function but with the
            // EquationSystems::READ_LEGACY_FORMAT bit set.
            this->read (name, mode, (read_flags | EquationSystems::READ_LEGACY_FORMAT));
            return;
          }

        // Figure out the libMesh version that created this file
        std::istringstream iss(version.substr(lm_pos + 8));
        int ver_major = 0, ver_minor = 0, ver_patch = 0;
        char dot;
        iss >> ver_major >> dot >> ver_minor >> dot >> ver_patch;
        io.set_version(LIBMESH_VERSION_ID(ver_major, ver_minor, ver_patch));


        read_parallel_files = (version.rfind(" parallel") < version.size());

        // If requested that we try to read infinite element information,
        // and the string " with infinite elements" is not in the version,
        // then tack it on.  This is for compatibility reading ifem
        // files written prior to 11/10/2008 - BSK
        if (try_read_ifems)
          if (!(version.rfind(" with infinite elements") < version.size()))
            version += " with infinite elements";

      }
    else
      libmesh_deprecated();

  START_LOG("read()","EquationSystems");

    // 2.)
    // Read the number of equation systems
    unsigned int n_sys=0;
    if (libMesh::processor_id() == 0) io.data (n_sys);
    CommWorld.broadcast(n_sys);

    for (unsigned int sys=0; sys<n_sys; sys++)
      {
        // 3.)
        // Read the name of the sys-th equation system
        std::string sys_name;
        if (libMesh::processor_id() == 0) io.data (sys_name);
        CommWorld.broadcast(sys_name);

        // 4.)
        // Read the type of the sys-th equation system
        std::string sys_type;
        if (libMesh::processor_id() == 0) io.data (sys_type);
        CommWorld.broadcast(sys_type);

        if (read_header)
          this->add_system (sys_type, sys_name);

        // 5.) - 9.)
        // Let System::read_header() do the job
        System& new_system = this->get_system(sys_name);
        new_system.read_header (io,
                                version,
                                read_header,
                                read_additional_data,
                                read_legacy_format);

        xda_systems.insert(std::make_pair(sys_name, &new_system));

        // If we're only creating "basic" systems, we need to tell
        // each system that before we call init() later.
        if (read_basic_only)
          new_system.set_basic_system_only();
      }
  }



  // Now we are ready to initialize the underlying data
  // structures. This will initialize the vectors for
  // storage, the dof_map, etc...
  if (read_header)
    this->init();

  // 10.) & 11.)
  // Read and set the numeric vector values
  if (read_data)
    {
      // the EquationSystems::read() method should look constant from the mesh
      // perspective, but we need to assign a temporary numbering to the nodes
      // and elements in the mesh, which requires that we abuse const_cast
      if (!read_legacy_format)
        {
          MeshBase &mesh = const_cast<MeshBase&>(this->get_mesh());
          MeshTools::Private::globally_renumber_nodes_and_elements(mesh);
        }

      Xdr local_io (read_parallel_files ? local_file_name(name) : "", mode);

      std::map<std::string, System*>::iterator
        pos = xda_systems.begin();

      for (; pos != xda_systems.end(); ++pos)
        if (read_legacy_format)
          {
            libmesh_deprecated();
            pos->second->read_legacy_data (io, read_additional_data);
          }
        else
          if (read_parallel_files)
            pos->second->read_parallel_data   (local_io, read_additional_data);
          else
            pos->second->read_serialized_data (io, read_additional_data);


      // Undo the temporary numbering.
      if (!read_legacy_format)
        _mesh.fix_broken_node_and_element_numbering();
    }

  STOP_LOG("read()","EquationSystems");

  // Localize each system's data
  this->update();
}



void EquationSystems::write(const std::string& name,
                            const unsigned int write_flags) const
{
  libMeshEnums::XdrMODE mode = WRITE;
  if (name.find(".xdr") != std::string::npos)
    mode = ENCODE;
  this->write(name, mode, write_flags);
}



void EquationSystems::write(const std::string& name,
                            const libMeshEnums::XdrMODE mode,
                            const unsigned int write_flags) const
{
  // the EquationSystems::write() method should look constant,
  // but we need to assign a temporary numbering to the nodes
  // and elements in the mesh, which requires that we abuse const_cast
  {
    MeshBase &mesh = const_cast<MeshBase&>(this->get_mesh());
    MeshTools::Private::globally_renumber_nodes_and_elements(mesh);
  }

   // set booleans from write_flags argument
   const bool write_data            = write_flags & EquationSystems::WRITE_DATA;
   const bool write_additional_data = write_flags & EquationSystems::WRITE_ADDITIONAL_DATA;

   // always write parallel files if we're instructed to write in
   // parallel
   const bool write_parallel_files  =
     (write_flags & EquationSystems::WRITE_PARALLEL_FILES) ||
   // but also write parallel files if we haven't been instructed to
   // write in serial and we're on a distributed mesh
     (!(write_flags & EquationSystems::WRITE_SERIAL_FILES) &&
      !this->get_mesh().is_serial());

  // New scope so that io will close before we try to zip the file
  {
    Xdr io((libMesh::processor_id()==0) ? name : "", mode);
    libmesh_assert (io.writing());

    START_LOG("write()","EquationSystems");

    const unsigned int proc_id = libMesh::processor_id();
    unsigned int n_sys         = this->n_systems();

    std::map<std::string, System*>::const_iterator
      pos = _systems.begin();

    std::string comment;
    char buf[256];

    // set the version number in the Xdr object
    io.set_version(LIBMESH_VERSION_ID(LIBMESH_MAJOR_VERSION,
                                      LIBMESH_MINOR_VERSION,
                                      LIBMESH_MICRO_VERSION));
    
    // Only write the header information
    // if we are processor 0.
    if (proc_id == 0)
      {
        // 1.)
        // Write the version header
        std::string version("libMesh-" + libMesh::get_io_compatibility_version());
        if (write_parallel_files) version += " parallel";

#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
        version += " with infinite elements";
#endif
        io.data (version, "# File Format Identifier");

        // 2.)
        // Write the number of equation systems
        io.data (n_sys, "# No. of Equation Systems");

        while (pos != _systems.end())
          {
            // 3.)
            // Write the name of the sys_num-th system
            {
              const unsigned int sys_num = pos->second->number();
              std::string sys_name       = pos->first;

              comment =  "# Name, System No. ";
              std::sprintf(buf, "%d", sys_num);
              comment += buf;

              io.data (sys_name, comment.c_str());
            }

            // 4.)
            // Write the type of system handled
            {
              const unsigned int sys_num = pos->second->number();
              std::string sys_type       = pos->second->system_type();

              comment =  "# Type, System No. ";
              std::sprintf(buf, "%d", sys_num);
              comment += buf;

              io.data (sys_type, comment.c_str());
            }

            // 5.) - 9.)
            // Let System::write_header() do the job
            pos->second->write_header (io, version, write_additional_data);

            ++pos;
          }
      }

    // Start from the first system, again,
    // to write vectors to disk, if wanted
    if (write_data)
      {
        // open a parallel buffer if warranted.
        Xdr local_io (write_parallel_files ? local_file_name(name) : "", mode);

        for (pos = _systems.begin(); pos != _systems.end(); ++pos)
          {
            // 10.) + 11.)
            if (write_parallel_files)
              pos->second->write_parallel_data (local_io,write_additional_data);
            else
              pos->second->write_serialized_data (io,write_additional_data);
          }
      }

    STOP_LOG("write()","EquationSystems");
  }

  // the EquationSystems::write() method should look constant,
  // but we need to undo the temporary numbering of the nodes
  // and elements in the mesh, which requires that we abuse const_cast
  const_cast<MeshBase&>(_mesh).fix_broken_node_and_element_numbering();
}

} // namespace libMesh

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC

Hosted By: