io_example.cpp

/*------------------------------------------------------------------------*/
/*                 Copyright 2010, 2011 Sandia Corporation.                     */
/*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
/*  license for use of this work by or on behalf of the U.S. Government.  */
/*  Export of this program may require a license from the                 */
/*  United States Government.                                             */
/*------------------------------------------------------------------------*/
 
#include <iostream>
#include <assert.h>
 
#include <stk_util/parallel/Parallel.hpp>
#include <init/Ionit_Initializer.h>
#include <Ioss_SubSystem.h>
 
#include <stk_mesh/base/Field.hpp>
#include <stk_mesh/base/FieldData.hpp>
#include <stk_mesh/base/BulkData.hpp>
 
#include <stk_mesh/fem/FEMMetaData.hpp>
#include <stk_mesh/fem/TopologyDimensions.hpp>
 
#include <stk_io/IossBridge.hpp>
 
namespace stk_example_io {
 
  void process_nodeblocks    (Ioss::Region &region, stk_classic::mesh::MetaData &meta);
 
  void process_elementblocks (Ioss::Region &region, stk_classic::mesh::MetaData &meta);
 
  void process_nodesets      (Ioss::Region &region, stk_classic::mesh::MetaData &meta);
 
  void process_sidesets      (Ioss::Region &region, stk_classic::mesh::MetaData &meta);
 
  void process_nodeblocks    (Ioss::Region &region, stk_classic::mesh::BulkData &bulk);
 
  void process_elementblocks (Ioss::Region &region, stk_classic::mesh::BulkData &bulk);
 
  void process_nodesets      (Ioss::Region &region, stk_classic::mesh::BulkData &bulk);
 
  void process_sidesets      (Ioss::Region &region, stk_classic::mesh::BulkData &bulk);
 
  void process_input_request (Ioss::Region &region, stk_classic::mesh::BulkData &bulk, int step);
 
  void process_output_request(Ioss::Region &region, stk_classic::mesh::BulkData &bulk, int step);
 
  void io_example( stk_classic::ParallelMachine comm,
       const std::string& in_filename,
       const std::string& out_filename,
       const std::string& decomp_method)
  {
    // Initialize IO system.  Registers all element types and storage
    // types and the exodusII default database type.
    Ioss::Init::Initializer init_db;
 
    std::cout << "========================================================================\n"
        << " Copy input mesh to output mesh.                                        \n"
        << "========================================================================\n";
 
    std::string dbtype("exodusII");
    Ioss::PropertyManager properties;
    if (!decomp_method.empty()) {
      properties.add(Ioss::Property("DECOMPOSITION_METHOD", Ioss::Utils::uppercase(decomp_method)));
    }
    Ioss::DatabaseIO *dbi = Ioss::IOFactory::create(dbtype, in_filename, Ioss::READ_MODEL,
                comm, properties);
    if (dbi == NULL || !dbi->ok()) {
      std::cerr  << "ERROR: Could not open database '" << in_filename
     << "' of type '" << dbtype << "'\n";
      std::exit(EXIT_FAILURE);
    }
 
    std::cout << "Reading input file:   " << in_filename << "\n";
    // NOTE: 'in_region' owns 'dbi' pointer at this time...
    Ioss::Region in_region(dbi, "input_model");
 
    // SUBSETTING PARSING/PREPROCESSING...
    // Just an example of how application could control whether an
    // entity is subsetted or not...
 
 
#if 0
    // Example command line in current code corresponding to behavior below:
    std::cout << "\nWhen processing file multi-block.g for use case 2, the blocks below will be omitted:\n";
    std::cout << "\tOMIT BLOCK Cblock Eblock I1 I2\n\n";
    Ioss::ElementBlock *eb = in_region.get_element_block("cblock");
    if (eb != NULL)
      eb->property_add(Ioss::Property(std::string("omitted"), 1));
 
    eb = in_region.get_element_block("eblock");
    if (eb != NULL)
      eb->property_add(Ioss::Property(std::string("omitted"), 1));
 
    eb = in_region.get_element_block("i1");
    if (eb != NULL)
      eb->property_add(Ioss::Property(std::string("omitted"), 1));
 
    eb = in_region.get_element_block("i2");
    if (eb != NULL)
      eb->property_add(Ioss::Property(std::string("omitted"), 1));
#endif
 
#if 0
    // Example for subsetting -- omit "odd" blocks
    if (entity->type() == Ioss::ELEMENTBLOCK) {
      int id = entity->get_property("id").get_int();
      if (id % 2) {
  entity->property_add(Ioss::Property(std::string("omitted"), 1));
  std::cout << "Skipping " << entity->type_string() << ": "  << entity->name() << "\n";
      }
    }
#endif
 
    //----------------------------------
    // Process Entity Types. Subsetting is possible.
 
    static size_t spatial_dimension = in_region.get_property("spatial_dimension").get_int();
 
    stk_classic::mesh::fem::FEMMetaData fem_meta_data( spatial_dimension );
    stk_classic::mesh::MetaData &meta_data = fem_meta_data.get_meta_data(fem_meta_data);
    process_elementblocks(in_region, meta_data);
    process_nodeblocks(in_region,    meta_data);
    process_sidesets(in_region,      meta_data);
    process_nodesets(in_region,      meta_data);
 
    //----------------------------------
    // Done populating meta data, commit and create bulk data
    meta_data.commit();
 
    //----------------------------------
    // Process Bulkdata for all Entity Types. Subsetting is possible.
    stk_classic::mesh::BulkData bulk_data(meta_data, comm);
 
    bulk_data.modification_begin();
    process_elementblocks(in_region, bulk_data);
    process_nodeblocks(in_region,    bulk_data);
    process_sidesets(in_region,      bulk_data);
    process_nodesets(in_region,      bulk_data);
    bulk_data.modification_end();
 
    //----------------------------------
    // OUTPUT...Create the output "mesh" portion
 
    std::cout << "Creating output file: " << out_filename << "\n";
    Ioss::DatabaseIO *dbo = Ioss::IOFactory::create(dbtype, out_filename,
                Ioss::WRITE_RESULTS,
                comm);
    if (dbo == NULL || !dbo->ok()) {
      std::cerr << "ERROR: Could not open results database '" << out_filename
    << "' of type '" << dbtype << "'\n";
      std::exit(EXIT_FAILURE);
    }
 
#if 0
    {
      // Code to test the remove_io_part_attribute functionality.
      // Hook this up to a command line option at some point to test nightly...
      const stk_classic::mesh::PartVector & all_parts = meta_data.get_parts();
      for ( stk_classic::mesh::PartVector::const_iterator ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
  stk_classic::mesh::Part * const part = *ip;
  const stk_classic::mesh::EntityRank part_rank = part->primary_entity_rank();
  
  if (stk_classic::io::is_part_io_part(*part) && part_rank == 2) {
    std::cout << "Removing part attribute from " << part->name() << "\n";
    stk_classic::io::remove_io_part_attribute(*part);
  }
      }
    }
#endif
 
    // NOTE: 'out_region' owns 'dbo' pointer at this time...
    Ioss::Region out_region(dbo, "results_output");
 
    stk_classic::io::define_output_db(out_region, bulk_data, &in_region);
    stk_classic::io::write_output_db(out_region,  bulk_data);
 
    // ------------------------------------------------------------------------
    out_region.begin_mode(Ioss::STATE_DEFINE_TRANSIENT);
 
    // Special processing for nodeblock (all nodes in model)...
    stk_classic::io::ioss_add_fields(meta_data.universal_part(), stk_classic::mesh::fem::FEMMetaData::NODE_RANK,
           out_region.get_node_blocks()[0],
           Ioss::Field::TRANSIENT);
 
    const stk_classic::mesh::PartVector & all_parts = meta_data.get_parts();
    for ( stk_classic::mesh::PartVector::const_iterator
      ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
 
      stk_classic::mesh::Part * const part = *ip;
 
      const stk_classic::mesh::EntityRank part_rank = part->primary_entity_rank();
 
      // Check whether this part should be output to results database.
      if (stk_classic::io::is_part_io_part(*part)) {
  // Get Ioss::GroupingEntity corresponding to this part...
  Ioss::GroupingEntity *entity = out_region.get_entity(part->name());
  if (entity != NULL) {
    if (entity->type() == Ioss::SIDESET) {
      Ioss::SideSet *sset = dynamic_cast<Ioss::SideSet*>(entity);
      assert(sset != NULL);
      int block_count = sset->block_count();
      for (int i=0; i < block_count; i++) {
        Ioss::SideBlock *fb = sset->get_block(i);
        stk_classic::io::ioss_add_fields(*part, part_rank,
               fb, Ioss::Field::TRANSIENT);
      }
    } else {
      stk_classic::io::ioss_add_fields(*part, part_rank,
             entity, Ioss::Field::TRANSIENT);
    }
  } else {
  }
      }
    }
    out_region.end_mode(Ioss::STATE_DEFINE_TRANSIENT);
    // ------------------------------------------------------------------------
 
    // Read and Write transient fields...
    out_region.begin_mode(Ioss::STATE_TRANSIENT);
    int timestep_count = in_region.get_property("state_count").get_int();
    for (int step = 1; step <= timestep_count; step++) {
      double time = in_region.get_state_time(step);
 
      // Read data from the io input mesh database into stk_classic::mesh fields...
      process_input_request(in_region, bulk_data, step);
 
      // execute()
 
      // Write data from the stk_classic::mesh fields out to the output database.a
      int out_step = out_region.add_state(time);
      process_output_request(out_region, bulk_data, out_step);
    }
    out_region.end_mode(Ioss::STATE_TRANSIENT);
  }
 
  // ========================================================================
  void process_nodeblocks(Ioss::Region &region, stk_classic::mesh::MetaData &meta)
  {
    const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
    assert(node_blocks.size() == 1);
 
    Ioss::NodeBlock *nb = node_blocks[0];
 
    assert(nb->field_exists("mesh_model_coordinates"));
    Ioss::Field coordinates = nb->get_field("mesh_model_coordinates");
    int spatial_dim = coordinates.transformed_storage()->component_count();
 
    stk_classic::mesh::Field<double,stk_classic::mesh::Cartesian> & coord_field =
      meta.declare_field<stk_classic::mesh::Field<double,stk_classic::mesh::Cartesian> >("coordinates");
 
    stk_classic::mesh::put_field( coord_field, stk_classic::mesh::fem::FEMMetaData::NODE_RANK, meta.universal_part(),
                          spatial_dim);
 
    stk_classic::io::define_io_fields(nb, Ioss::Field::TRANSIENT, meta.universal_part(),stk_classic::mesh::fem::FEMMetaData::NODE_RANK);
  }
 
  // ========================================================================
  void process_elementblocks(Ioss::Region &region, stk_classic::mesh::MetaData &meta)
  {
    const stk_classic::mesh::fem::FEMMetaData &fem_meta_data = stk_classic::mesh::fem::FEMMetaData::get(meta);
    const stk_classic::mesh::EntityRank element_rank = fem_meta_data.element_rank();
 
    const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
    stk_classic::io::default_part_processing(elem_blocks, meta, element_rank);
 
    // Parts were created above, now handle element block specific
    // information (topology, attributes, ...);
    for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
  it != elem_blocks.end(); ++it) {
      Ioss::ElementBlock *entity = *it;
 
      if (stk_classic::io::include_entity(entity)) {
  stk_classic::mesh::Part* const part = meta.get_part(entity->name());
  assert(part != NULL);
 
      const stk_classic::mesh::EntityRank part_rank = part->primary_entity_rank();
 
  // Element Block attributes (if any)...
  stk_classic::io::define_io_fields(entity, Ioss::Field::ATTRIBUTE,
          *part,
                                  part_rank);
 
  stk_classic::io::define_io_fields(entity, Ioss::Field::TRANSIENT,
          *part,
                                  part_rank);
 
  const CellTopologyData* cell_topo = fem_meta_data.get_cell_topology(*part).getCellTopologyData();
  std::string cell_topo_name = "UNKNOWN";
  if (cell_topo != NULL)
    cell_topo_name = cell_topo->name;
      }
    }
  }
 
  // ========================================================================
  void process_nodesets(Ioss::Region &region, stk_classic::mesh::MetaData &meta)
  {
    const Ioss::NodeSetContainer& node_sets = region.get_nodesets();
    stk_classic::io::default_part_processing(node_sets, meta, stk_classic::mesh::fem::FEMMetaData::NODE_RANK);
 
    stk_classic::mesh::Field<double> & distribution_factors_field =
      meta.declare_field<stk_classic::mesh::Field<double> >("distribution_factors");
 
    for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
  it != node_sets.end(); ++it) {
      Ioss::NodeSet *entity = *it;
 
      if (stk_classic::io::include_entity(entity)) {
  stk_classic::mesh::Part* const part = meta.get_part(entity->name());
  assert(part != NULL);
  assert(entity->field_exists("distribution_factors"));
 
  stk_classic::mesh::put_field(distribution_factors_field, stk_classic::mesh::fem::FEMMetaData::NODE_RANK, *part);
 
  stk_classic::io::define_io_fields(entity, Ioss::Field::TRANSIENT,
          *part,
                                  part->primary_entity_rank() ) ;
      }
    }
  }
 
  // ========================================================================
  void process_surface_entity(Ioss::SideSet *sset, stk_classic::mesh::MetaData &meta,
            stk_classic::mesh::EntityRank sset_rank)
  {
    assert(sset->type() == Ioss::SIDESET);
    Ioss::SideSet *fs = dynamic_cast<Ioss::SideSet *>(sset);
    assert(fs != NULL);
    const Ioss::SideBlockContainer& blocks = fs->get_side_blocks();
    stk_classic::io::default_part_processing(blocks, meta, sset_rank);
 
    stk_classic::mesh::Part* const fs_part = meta.get_part(sset->name());
    assert(fs_part != NULL);
 
    stk_classic::mesh::Field<double, stk_classic::mesh::ElementNode> *distribution_factors_field = NULL;
    bool surface_df_defined = false; // Has the surface df field been defined yet?
 
 
    int block_count = sset->block_count();
    for (int i=0; i < block_count; i++) {
      Ioss::SideBlock *side_block = sset->get_block(i);
      if (stk_classic::io::include_entity(side_block)) {
  stk_classic::mesh::Part * const side_block_part = meta.get_part(side_block->name());
  assert(side_block_part != NULL);
  meta.declare_part_subset(*fs_part, *side_block_part);
 
        const stk_classic::mesh::EntityRank part_rank = side_block_part->primary_entity_rank();
 
  if (side_block->field_exists("distribution_factors")) {
    if (!surface_df_defined) {
      std::string field_name = sset->name() + "_distribution_factors";
      distribution_factors_field =
        &meta.declare_field<stk_classic::mesh::Field<double, stk_classic::mesh::ElementNode> >(field_name);
      stk_classic::io::set_distribution_factor_field(*fs_part, *distribution_factors_field);
      surface_df_defined = true;
    }
    stk_classic::io::set_distribution_factor_field(*side_block_part, *distribution_factors_field);
    int side_node_count = side_block->topology()->number_nodes();
    stk_classic::mesh::put_field(*distribution_factors_field,
                               part_rank,
                               *side_block_part, side_node_count);
  }
 
  stk_classic::io::define_io_fields(side_block, Ioss::Field::TRANSIENT,
          *side_block_part,
                                  part_rank);
      }
    }
  }
 
  // ========================================================================
  void process_sidesets(Ioss::Region &region, stk_classic::mesh::MetaData &meta)
  {
    stk_classic::mesh::fem::FEMMetaData &fem = stk_classic::mesh::fem::FEMMetaData::get(meta);
    const stk_classic::mesh::EntityRank side_rank = fem.side_rank();
 
    const Ioss::SideSetContainer& side_sets = region.get_sidesets();
    stk_classic::io::default_part_processing(side_sets, meta, side_rank);
 
    for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
  it != side_sets.end(); ++it) {
      Ioss::SideSet *entity = *it;
 
      if (stk_classic::io::include_entity(entity)) {
  process_surface_entity(entity, meta, side_rank);
      }
    }
  }
 
  // ========================================================================
  // Bulk Data
  // ========================================================================
  void process_nodeblocks(Ioss::Region &region, stk_classic::mesh::BulkData &bulk)
  {
    // This must be called after the "process_element_blocks" call
    // since there may be nodes that exist in the database that are
    // not part of the analysis mesh due to subsetting of the element
    // blocks.
 
    const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
    assert(node_blocks.size() == 1);
 
    Ioss::NodeBlock *nb = node_blocks[0];
 
    std::vector<stk_classic::mesh::Entity*> nodes;
    stk_classic::io::get_entity_list(nb, stk_classic::mesh::fem::FEMMetaData::NODE_RANK, bulk, nodes);
 
    const stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(bulk);
    stk_classic::mesh::Field<double,stk_classic::mesh::Cartesian> *coord_field =
      meta.get_field<stk_classic::mesh::Field<double,stk_classic::mesh::Cartesian> >("coordinates");
 
    stk_classic::io::field_data_from_ioss(coord_field, nodes, nb, "mesh_model_coordinates");
  }
 
  // ========================================================================
  void process_elementblocks(Ioss::Region &region, stk_classic::mesh::BulkData &bulk)
  {
    const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
 
    for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
  it != elem_blocks.end(); ++it) {
      Ioss::ElementBlock *entity = *it;
 
      if (stk_classic::io::include_entity(entity)) {
  const std::string &name = entity->name();
        const stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(bulk);
        const stk_classic::mesh::fem::FEMMetaData &fem_meta_data = stk_classic::mesh::fem::FEMMetaData::get(meta);
  stk_classic::mesh::Part* const part = meta.get_part(name);
  assert(part != NULL);
 
  const CellTopologyData* cell_topo = fem_meta_data.get_cell_topology(*part).getCellTopologyData();
  if (cell_topo == NULL) {
          std::ostringstream msg ;
    msg << " INTERNAL_ERROR: Part " << part->name() << " returned NULL from get_cell_topology()";
    throw std::runtime_error( msg.str() );
  }
 
  std::vector<int> elem_ids ;
  std::vector<int> connectivity ;
  std::vector<stk_classic::mesh::EntityId> connectivity2 ;
 
  entity->get_field_data("ids", elem_ids);
  entity->get_field_data("connectivity", connectivity);
        connectivity2.reserve(connectivity.size());
        std::copy(connectivity.begin(), connectivity.end(), std::back_inserter(connectivity2));
 
  size_t element_count = elem_ids.size();
  int nodes_per_elem = cell_topo->node_count ;
 
  std::vector<stk_classic::mesh::Entity*> elements(element_count);
  for(size_t i=0; i<element_count; ++i) {
    stk_classic::mesh::EntityId *conn = &connectivity2[i*nodes_per_elem];
    elements[i] = &stk_classic::mesh::fem::declare_element(bulk, *part, elem_ids[i], conn);
  }
 
  // For this example, we are just taking all attribute fields
  // found on the io database and populating fields on the
  // corresponding mesh part.  In practice, would probably be
  // selective about which attributes to use...
  Ioss::NameList names;
  entity->field_describe(Ioss::Field::ATTRIBUTE, &names);
  for (Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
    if (*I == "attribute" && names.size() > 1)
      continue;
    stk_classic::mesh::FieldBase *field = meta.get_field<stk_classic::mesh::FieldBase>(*I);
    stk_classic::io::field_data_from_ioss(field, elements, entity, *I);
 
  }
      }
    }
  }
 
  // ========================================================================
  void process_nodesets(Ioss::Region &region, stk_classic::mesh::BulkData &bulk)
  {
    // Should only process nodes that have already been defined via the element
    // blocks connectivity lists.
    const Ioss::NodeSetContainer& node_sets = region.get_nodesets();
 
    for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
  it != node_sets.end(); ++it) {
      Ioss::NodeSet *entity = *it;
 
      if (stk_classic::io::include_entity(entity)) {
  const std::string & name = entity->name();
  const stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(bulk);
  stk_classic::mesh::Part* const part = meta.get_part(name);
  assert(part != NULL);
  stk_classic::mesh::PartVector add_parts( 1 , part );
 
  std::vector<int> node_ids ;
  int node_count = entity->get_field_data("ids", node_ids);
 
  std::vector<stk_classic::mesh::Entity*> nodes(node_count);
  for(int i=0; i<node_count; ++i) {
    nodes[i] = bulk.get_entity( stk_classic::mesh::fem::FEMMetaData::NODE_RANK, node_ids[i] );
    if (nodes[i] != NULL)
      bulk.declare_entity(stk_classic::mesh::fem::FEMMetaData::NODE_RANK, node_ids[i], add_parts );
  }
 
  stk_classic::mesh::Field<double> *df_field =
    meta.get_field<stk_classic::mesh::Field<double> >("distribution_factors");
 
  if (df_field != NULL) {
    stk_classic::io::field_data_from_ioss(df_field, nodes, entity, "distribution_factors");
  }
      }
    }
  }
 
  // ========================================================================
  void process_surface_entity(const Ioss::SideSet* sset ,
            stk_classic::mesh::BulkData & bulk)
  {
    assert(sset->type() == Ioss::SIDESET);
 
    const stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(bulk);
    const stk_classic::mesh::fem::FEMMetaData &fem_meta_data = stk_classic::mesh::fem::FEMMetaData::get(meta);
    const stk_classic::mesh::EntityRank element_rank = fem_meta_data.element_rank();
 
    int block_count = sset->block_count();
    for (int i=0; i < block_count; i++) {
      Ioss::SideBlock *block = sset->get_block(i);
      if (stk_classic::io::include_entity(block)) {
  std::vector<int> side_ids ;
  std::vector<int> elem_side ;
 
  stk_classic::mesh::Part * const side_block_part = meta.get_part(block->name());
  stk_classic::mesh::EntityRank side_rank = side_block_part->primary_entity_rank();
 
  block->get_field_data("ids", side_ids);
  block->get_field_data("element_side", elem_side);
 
  assert(side_ids.size() * 2 == elem_side.size());
  stk_classic::mesh::PartVector add_parts( 1 , side_block_part );
 
  size_t side_count = side_ids.size();
  std::vector<stk_classic::mesh::Entity*> sides(side_count);
  for(size_t is=0; is<side_count; ++is) {
 
    stk_classic::mesh::Entity* const elem = bulk.get_entity(element_rank, elem_side[is*2]);
 
    // If NULL, then the element was probably assigned to an
    // element block that appears in the database, but was
    // subsetted out of the analysis mesh. Only process if
    // non-null.
    if (elem != NULL) {
      // Ioss uses 1-based side ordinal, stk_classic::mesh uses 0-based.
      // Hence the '-1' in the following line.
      int side_ordinal = elem_side[is*2+1] - 1 ;
 
      stk_classic::mesh::Entity *side = NULL;
      if (side_rank == 2) {
        side = &stk_classic::mesh::fem::declare_element_side(bulk, side_ids[is], *elem, side_ordinal);
      } else {
        side = &stk_classic::mesh::fem::declare_element_edge(bulk, side_ids[is], *elem, side_ordinal);
      }
      bulk.change_entity_parts( *side, add_parts );
      sides[is] = side;
    } else {
      sides[is] = NULL;
    }
  }
 
  const stk_classic::mesh::Field<double, stk_classic::mesh::ElementNode> *df_field =
    stk_classic::io::get_distribution_factor_field(*side_block_part);
  if (df_field != NULL) {
    stk_classic::io::field_data_from_ioss(df_field, sides, block, "distribution_factors");
  }
      }
    }
  }
 
  // ========================================================================
  void process_sidesets(Ioss::Region &region, stk_classic::mesh::BulkData &bulk)
  {
    const Ioss::SideSetContainer& side_sets = region.get_sidesets();
 
    for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
  it != side_sets.end(); ++it) {
      Ioss::SideSet *entity = *it;
 
      if (stk_classic::io::include_entity(entity)) {
  process_surface_entity(entity, bulk);
      }
    }
  }
 
  // ========================================================================
  // ========================================================================
  void get_field_data(stk_classic::mesh::BulkData &bulk, stk_classic::mesh::Part &part,
          stk_classic::mesh::EntityRank part_type,
          Ioss::GroupingEntity *io_entity,
          Ioss::Field::RoleType filter_role)
  {
    std::vector<stk_classic::mesh::Entity*> entities;
    stk_classic::io::get_entity_list(io_entity, part_type, bulk, entities);
 
    stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(part);
    stk_classic::mesh::Part &universal = meta.universal_part();
    const std::vector<stk_classic::mesh::FieldBase*> &fields = meta.get_fields();
 
    std::vector<stk_classic::mesh::FieldBase *>::const_iterator I = fields.begin();
    while (I != fields.end()) {
      const stk_classic::mesh::FieldBase *f = *I; ++I;
      if (stk_classic::io::is_valid_part_field(f, part_type, part, universal, filter_role)) {
  stk_classic::io::field_data_from_ioss(f, entities, io_entity, f->name());
      }
    }
  }
 
  void process_input_request(Ioss::Region &region,
           stk_classic::mesh::BulkData &bulk,
           int step)
  {
    region.begin_state(step);
 
    // Special processing for nodeblock (all nodes in model)...
    const stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(bulk);
 
    // ??? Get field data from nodeblock...
    get_field_data(bulk, meta.universal_part(), stk_classic::mesh::fem::FEMMetaData::NODE_RANK,
       region.get_node_blocks()[0], Ioss::Field::TRANSIENT);
 
    const stk_classic::mesh::PartVector & all_parts = meta.get_parts();
    for ( stk_classic::mesh::PartVector::const_iterator
      ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
 
      stk_classic::mesh::Part * const part = *ip;
 
      const stk_classic::mesh::EntityRank part_rank = part->primary_entity_rank();
 
      // Check whether this part should be output to results database.
      if (stk_classic::io::is_part_io_part(*part)) {
  // Get Ioss::GroupingEntity corresponding to this part...
  Ioss::GroupingEntity *entity = region.get_entity(part->name());
  if (entity != NULL) {
    if (entity->type() == Ioss::SIDESET) {
      Ioss::SideSet *sset = dynamic_cast<Ioss::SideSet*>(entity);
      assert(sset != NULL);
      int block_count = sset->block_count();
      for (int i=0; i < block_count; i++) {
        Ioss::SideBlock *side_block = sset->get_block(i);
    get_field_data(bulk, *part,
                               part_rank,
             side_block, Ioss::Field::TRANSIENT);
      }
    } else {
      get_field_data(bulk, *part,
                           part_rank,
         entity, Ioss::Field::TRANSIENT);
    }
  } else {
  }
      }
    }
 
    region.end_state(step);
  }
 
  void put_field_data(stk_classic::mesh::BulkData &bulk, stk_classic::mesh::Part &part,
          stk_classic::mesh::EntityRank part_type,
          Ioss::GroupingEntity *io_entity,
          Ioss::Field::RoleType filter_role)
  {
    std::vector<stk_classic::mesh::Entity*> entities;
    stk_classic::io::get_entity_list(io_entity, part_type, bulk, entities);
 
    stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(part);
    stk_classic::mesh::Part &universal = meta.universal_part();
    const std::vector<stk_classic::mesh::FieldBase*> &fields = meta.get_fields();
 
    std::vector<stk_classic::mesh::FieldBase *>::const_iterator I = fields.begin();
    while (I != fields.end()) {
      const stk_classic::mesh::FieldBase *f = *I; ++I;
      if (stk_classic::io::is_valid_part_field(f, part_type, part, universal, filter_role)) {
  stk_classic::io::field_data_to_ioss(f, entities, io_entity, f->name(), filter_role);
      }
    }
  }
 
  void process_output_request(Ioss::Region &region,
            stk_classic::mesh::BulkData &bulk,
            int step)
  {
    region.begin_state(step);
    // Special processing for nodeblock (all nodes in model)...
    const stk_classic::mesh::MetaData& meta = stk_classic::mesh::MetaData::get(bulk);
 
    put_field_data(bulk, meta.universal_part(), stk_classic::mesh::fem::FEMMetaData::NODE_RANK,
       region.get_node_blocks()[0], Ioss::Field::TRANSIENT);
 
    const stk_classic::mesh::PartVector & all_parts = meta.get_parts();
    for ( stk_classic::mesh::PartVector::const_iterator
      ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
 
      stk_classic::mesh::Part * const part = *ip;
 
      const stk_classic::mesh::EntityRank part_rank = part->primary_entity_rank();
 
      // Check whether this part should be output to results database.
      if (stk_classic::io::is_part_io_part(*part)) {
 
  // Get Ioss::GroupingEntity corresponding to this part...
  Ioss::GroupingEntity *entity = region.get_entity(part->name());
  if (entity != NULL) {
 
    if (entity->type() == Ioss::SIDESET) {
      Ioss::SideSet *sset = dynamic_cast<Ioss::SideSet*>(entity);
      assert(sset != NULL);
      int block_count = sset->block_count();
 
      for (int i=0; i < block_count; i++) {
        Ioss::SideBlock *side_block = sset->get_block(i);
    put_field_data(bulk, *part, part_rank,
             side_block, Ioss::Field::TRANSIENT);
      }
    } else {
      put_field_data(bulk, *part, part_rank,
         entity, Ioss::Field::TRANSIENT);
    }
  } else {
  }
      }
    }
    region.end_state(step);
  }
}
 
// ========================================================================
#include <boost/program_options.hpp>
 
#include <stk_util/parallel/BroadcastArg.hpp>
#include <stk_util/environment/ProgramOptions.hpp>
 
  namespace bopt = boost::program_options;
  int main(int argc, char** argv)
  {
    //----------------------------------
    // Broadcast argc and argv to all processors.
 
    stk_classic::ParallelMachine comm = stk_classic::parallel_machine_init(&argc, &argv);
 
    stk_classic::BroadcastArg b_arg(comm, argc, argv);
 
    //----------------------------------
    // Process the broadcast command line arguments
 
    bopt::options_description desc("options");
 
    desc.add_options()
      ("help,h",        "produce help message")
      ("mesh",         bopt::value<std::string>(), "mesh file" )
      ("decomposition,D", bopt::value<std::string>(), "decomposition method" )
      ("directory,d",  bopt::value<std::string>(), "working directory" )
      ("output-log,o", bopt::value<std::string>(), "output log path" )
      ("runtest,r",    bopt::value<std::string>(), "runtest pid file" );
 
    stk_classic::get_options_description().add(desc);
 
    bopt::variables_map &vm = stk_classic::get_variables_map();
    try {
      bopt::store(bopt::parse_command_line(b_arg.m_argc, b_arg.m_argv, desc), vm);
      bopt::notify(vm);
    }
    catch (std::exception & /* x */) {
      std::exit(1);
    }
 
    if (vm.count("help")) {
      std::cout << desc << "\n";
      std::exit(EXIT_SUCCESS);
    }
 
    //----------------------------------
 
    if ( vm.count("mesh") ) {
      std::string in_filename = boost::any_cast<std::string>(vm["mesh"].value());
      std::string out_filename = in_filename + ".out";
      std::string decomp_method;
      if (vm.count("decomposition")) {
  decomp_method = boost::any_cast<std::string>(vm["decomposition"].value());
      }
      stk_example_io::io_example(comm, in_filename, out_filename, decomp_method );
    } else {
      std::cout << "OPTION ERROR: The '--mesh <filename>' option is required!\n";
      std::exit(EXIT_FAILURE);
    }
    stk_classic::parallel_machine_finalize();
 
    return 0;
  }