Sierra Toolkit Version of the Day
io_example.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010, 2011 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 #include <iostream>
00010 #include <assert.h>
00011 
00012 #include <stk_util/parallel/Parallel.hpp>
00013 #include <init/Ionit_Initializer.h>
00014 #include <Ioss_SubSystem.h>
00015 
00016 #include <stk_mesh/base/Field.hpp>
00017 #include <stk_mesh/base/FieldData.hpp>
00018 #include <stk_mesh/base/BulkData.hpp>
00019 
00020 #include <stk_mesh/fem/FEMMetaData.hpp>
00021 #include <stk_mesh/fem/TopologyDimensions.hpp>
00022 
00023 #include <stk_io/IossBridge.hpp>
00024 
00041 namespace stk_example_io {
00042 
00046   void process_nodeblocks    (Ioss::Region &region, stk::mesh::MetaData &meta);
00047 
00053   void process_elementblocks (Ioss::Region &region, stk::mesh::MetaData &meta);
00054 
00061   void process_nodesets      (Ioss::Region &region, stk::mesh::MetaData &meta);
00062 
00076   void process_sidesets      (Ioss::Region &region, stk::mesh::MetaData &meta);
00077 
00084   void process_nodeblocks    (Ioss::Region &region, stk::mesh::BulkData &bulk);
00085 
00097   void process_elementblocks (Ioss::Region &region, stk::mesh::BulkData &bulk);
00098 
00105   void process_nodesets      (Ioss::Region &region, stk::mesh::BulkData &bulk);
00106 
00112   void process_sidesets      (Ioss::Region &region, stk::mesh::BulkData &bulk);
00113 
00120   void process_input_request (Ioss::Region &region, stk::mesh::BulkData &bulk, int step);
00121 
00134   void process_output_request(Ioss::Region &region, stk::mesh::BulkData &bulk, int step);
00135 
00162   void io_example( stk::ParallelMachine comm,
00163        const std::string& in_filename,
00164        const std::string& out_filename)
00165   {
00166     // Initialize IO system.  Registers all element types and storage
00167     // types and the exodusII default database type.
00168     Ioss::Init::Initializer init_db;
00169 
00170     std::cout << "========================================================================\n"
00171         << " Copy input mesh to output mesh.                                        \n"
00172         << "========================================================================\n";
00173 
00174     std::string dbtype("exodusII");
00175     Ioss::DatabaseIO *dbi = Ioss::IOFactory::create(dbtype, in_filename, Ioss::READ_MODEL,
00176                 comm);
00177     if (dbi == NULL || !dbi->ok()) {
00178       std::cerr  << "ERROR: Could not open database '" << in_filename
00179      << "' of type '" << dbtype << "'\n";
00180       std::exit(EXIT_FAILURE);
00181     }
00182 
00183     std::cout << "Reading input file:   " << in_filename << "\n";
00184     // NOTE: 'in_region' owns 'dbi' pointer at this time...
00185     Ioss::Region in_region(dbi, "input_model");
00186 
00187     // SUBSETTING PARSING/PREPROCESSING...
00188     // Just an example of how application could control whether an
00189     // entity is subsetted or not...
00190 
00191 
00192 #if 0
00193     // Example command line in current code corresponding to behavior below:
00194     std::cout << "\nWhen processing file multi-block.g for use case 2, the blocks below will be omitted:\n";
00195     std::cout << "\tOMIT BLOCK Cblock Eblock I1 I2\n\n";
00196     Ioss::ElementBlock *eb = in_region.get_element_block("cblock");
00197     if (eb != NULL)
00198       eb->property_add(Ioss::Property(std::string("omitted"), 1));
00199 
00200     eb = in_region.get_element_block("eblock");
00201     if (eb != NULL)
00202       eb->property_add(Ioss::Property(std::string("omitted"), 1));
00203 
00204     eb = in_region.get_element_block("i1");
00205     if (eb != NULL)
00206       eb->property_add(Ioss::Property(std::string("omitted"), 1));
00207 
00208     eb = in_region.get_element_block("i2");
00209     if (eb != NULL)
00210       eb->property_add(Ioss::Property(std::string("omitted"), 1));
00211 #endif
00212 
00213 #if 0
00214     // Example for subsetting -- omit "odd" blocks
00215     if (entity->type() == Ioss::ELEMENTBLOCK) {
00216       int id = entity->get_property("id").get_int();
00217       if (id % 2) {
00218   entity->property_add(Ioss::Property(std::string("omitted"), 1));
00219   std::cout << "Skipping " << entity->type_string() << ": "  << entity->name() << "\n";
00220       }
00221     }
00222 #endif
00223 
00224     //----------------------------------
00225     // Process Entity Types. Subsetting is possible.
00226 
00227     static size_t spatial_dimension = in_region.get_property("spatial_dimension").get_int();
00228 
00229     stk::mesh::fem::FEMMetaData fem_meta_data( spatial_dimension );
00230     stk::mesh::MetaData &meta_data = fem_meta_data.get_meta_data(fem_meta_data);
00231     process_elementblocks(in_region, meta_data);
00232     process_nodeblocks(in_region,    meta_data);
00233     process_sidesets(in_region,      meta_data);
00234     process_nodesets(in_region,      meta_data);
00235 
00236     //----------------------------------
00237     // Done populating meta data, commit and create bulk data
00238     meta_data.commit();
00239 
00240     //----------------------------------
00241     // Process Bulkdata for all Entity Types. Subsetting is possible.
00242     stk::mesh::BulkData bulk_data(meta_data, comm);
00243 
00244     bulk_data.modification_begin();
00245     process_elementblocks(in_region, bulk_data);
00246     process_nodeblocks(in_region,    bulk_data);
00247     process_sidesets(in_region,      bulk_data);
00248     process_nodesets(in_region,      bulk_data);
00249     bulk_data.modification_end();
00250 
00251     //----------------------------------
00252     // OUTPUT...Create the output "mesh" portion
00253 
00254     std::cout << "Creating output file: " << out_filename << "\n";
00255     Ioss::DatabaseIO *dbo = Ioss::IOFactory::create(dbtype, out_filename,
00256                 Ioss::WRITE_RESULTS,
00257                 comm);
00258     if (dbo == NULL || !dbo->ok()) {
00259       std::cerr << "ERROR: Could not open results database '" << out_filename
00260     << "' of type '" << dbtype << "'\n";
00261       std::exit(EXIT_FAILURE);
00262     }
00263 
00264 #if 0
00265     {
00266       // Code to test the remove_io_part_attribute functionality.
00267       // Hook this up to a command line option at some point to test nightly...
00268       const stk::mesh::PartVector & all_parts = meta_data.get_parts();
00269       for ( stk::mesh::PartVector::const_iterator ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
00270   stk::mesh::Part * const part = *ip;
00271   const stk::mesh::EntityRank part_rank = part->primary_entity_rank();
00272   
00273   if (stk::io::is_part_io_part(*part) && part_rank == 2) {
00274     std::cout << "Removing part attribute from " << part->name() << "\n";
00275     stk::io::remove_io_part_attribute(*part);
00276   }
00277       }
00278     }
00279 #endif
00280 
00281     // NOTE: 'out_region' owns 'dbo' pointer at this time...
00282     Ioss::Region out_region(dbo, "results_output");
00283 
00284     stk::io::define_output_db(out_region, bulk_data, &in_region);
00285     stk::io::write_output_db(out_region,  bulk_data);
00286 
00287     // ------------------------------------------------------------------------
00302     out_region.begin_mode(Ioss::STATE_DEFINE_TRANSIENT);
00303 
00304     // Special processing for nodeblock (all nodes in model)...
00305     stk::io::ioss_add_fields(meta_data.universal_part(), stk::mesh::fem::FEMMetaData::NODE_RANK,
00306            out_region.get_node_blocks()[0],
00307            Ioss::Field::TRANSIENT);
00308 
00309     const stk::mesh::PartVector & all_parts = meta_data.get_parts();
00310     for ( stk::mesh::PartVector::const_iterator
00311       ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
00312 
00313       stk::mesh::Part * const part = *ip;
00314 
00315       const stk::mesh::EntityRank part_rank = part->primary_entity_rank();
00316 
00317       // Check whether this part should be output to results database.
00318       if (stk::io::is_part_io_part(*part)) {
00319   // Get Ioss::GroupingEntity corresponding to this part...
00320   Ioss::GroupingEntity *entity = out_region.get_entity(part->name());
00321   if (entity != NULL) {
00322     if (entity->type() == Ioss::SIDESET) {
00323       Ioss::SideSet *sset = dynamic_cast<Ioss::SideSet*>(entity);
00324       assert(sset != NULL);
00325       int block_count = sset->block_count();
00326       for (int i=0; i < block_count; i++) {
00327         Ioss::SideBlock *fb = sset->get_block(i);
00328         stk::io::ioss_add_fields(*part, part_rank,
00329                fb, Ioss::Field::TRANSIENT);
00330       }
00331     } else {
00332       stk::io::ioss_add_fields(*part, part_rank,
00333              entity, Ioss::Field::TRANSIENT);
00334     }
00335   } else {
00338   }
00339       }
00340     }
00341     out_region.end_mode(Ioss::STATE_DEFINE_TRANSIENT);
00342     // ------------------------------------------------------------------------
00343 
00344     // Read and Write transient fields...
00345     out_region.begin_mode(Ioss::STATE_TRANSIENT);
00346     int timestep_count = in_region.get_property("state_count").get_int();
00347     for (int step = 1; step <= timestep_count; step++) {
00348       double time = in_region.get_state_time(step);
00349 
00350       // Read data from the io input mesh database into stk::mesh fields...
00351       process_input_request(in_region, bulk_data, step);
00352 
00353       // execute()
00354 
00355       // Write data from the stk::mesh fields out to the output database.a
00356       int out_step = out_region.add_state(time);
00357       process_output_request(out_region, bulk_data, out_step);
00358     }
00359     out_region.end_mode(Ioss::STATE_TRANSIENT);
00360   }
00361 
00362   // ========================================================================
00363   void process_nodeblocks(Ioss::Region &region, stk::mesh::MetaData &meta)
00364   {
00365     const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
00366     assert(node_blocks.size() == 1);
00367 
00368     Ioss::NodeBlock *nb = node_blocks[0];
00369 
00370     assert(nb->field_exists("mesh_model_coordinates"));
00371     Ioss::Field coordinates = nb->get_field("mesh_model_coordinates");
00372     int spatial_dim = coordinates.transformed_storage()->component_count();
00373 
00374     stk::mesh::Field<double,stk::mesh::Cartesian> & coord_field =
00375       meta.declare_field<stk::mesh::Field<double,stk::mesh::Cartesian> >("coordinates");
00376 
00377     stk::mesh::put_field( coord_field, stk::mesh::fem::FEMMetaData::NODE_RANK, meta.universal_part(),
00378                           spatial_dim);
00379 
00384     stk::io::define_io_fields(nb, Ioss::Field::TRANSIENT, meta.universal_part(),stk::mesh::fem::FEMMetaData::NODE_RANK);
00385   }
00386 
00387   // ========================================================================
00388   void process_elementblocks(Ioss::Region &region, stk::mesh::MetaData &meta)
00389   {
00390     const stk::mesh::fem::FEMMetaData &fem_meta_data = stk::mesh::fem::FEMMetaData::get(meta);
00391     const stk::mesh::EntityRank element_rank = fem_meta_data.element_rank();
00392 
00393     const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
00394     stk::io::default_part_processing(elem_blocks, meta, element_rank);
00395 
00396     // Parts were created above, now handle element block specific
00397     // information (topology, attributes, ...);
00398     for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
00399   it != elem_blocks.end(); ++it) {
00400       Ioss::ElementBlock *entity = *it;
00401 
00402       if (stk::io::include_entity(entity)) {
00403   stk::mesh::Part* const part = meta.get_part(entity->name());
00404   assert(part != NULL);
00405 
00406       const stk::mesh::EntityRank part_rank = part->primary_entity_rank();
00407 
00408   // Element Block attributes (if any)...
00413   stk::io::define_io_fields(entity, Ioss::Field::ATTRIBUTE,
00414           *part,
00415                                   part_rank);
00416 
00421   stk::io::define_io_fields(entity, Ioss::Field::TRANSIENT,
00422           *part,
00423                                   part_rank);
00424 
00425   const CellTopologyData* cell_topo = fem_meta_data.get_cell_topology(*part).getCellTopologyData();
00426   std::string cell_topo_name = "UNKNOWN";
00427   if (cell_topo != NULL)
00428     cell_topo_name = cell_topo->name;
00429       }
00430     }
00431   }
00432 
00433   // ========================================================================
00434   void process_nodesets(Ioss::Region &region, stk::mesh::MetaData &meta)
00435   {
00436     const Ioss::NodeSetContainer& node_sets = region.get_nodesets();
00437     stk::io::default_part_processing(node_sets, meta, stk::mesh::fem::FEMMetaData::NODE_RANK);
00438 
00443     stk::mesh::Field<double> & distribution_factors_field =
00444       meta.declare_field<stk::mesh::Field<double> >("distribution_factors");
00445 
00451     for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
00452   it != node_sets.end(); ++it) {
00453       Ioss::NodeSet *entity = *it;
00454 
00455       if (stk::io::include_entity(entity)) {
00456   stk::mesh::Part* const part = meta.get_part(entity->name());
00457   assert(part != NULL);
00458   assert(entity->field_exists("distribution_factors"));
00459 
00460   stk::mesh::put_field(distribution_factors_field, stk::mesh::fem::FEMMetaData::NODE_RANK, *part);
00461 
00466   stk::io::define_io_fields(entity, Ioss::Field::TRANSIENT,
00467           *part,
00468                                   part->primary_entity_rank() ) ;
00469       }
00470     }
00471   }
00472 
00473   // ========================================================================
00474   void process_surface_entity(Ioss::SideSet *sset, stk::mesh::MetaData &meta,
00475             stk::mesh::EntityRank sset_rank)
00476   {
00477     assert(sset->type() == Ioss::SIDESET);
00478     Ioss::SideSet *fs = dynamic_cast<Ioss::SideSet *>(sset);
00479     assert(fs != NULL);
00480     const Ioss::SideBlockContainer& blocks = fs->get_side_blocks();
00481     stk::io::default_part_processing(blocks, meta, sset_rank);
00482 
00483     stk::mesh::Part* const fs_part = meta.get_part(sset->name());
00484     assert(fs_part != NULL);
00485 
00486     stk::mesh::Field<double, stk::mesh::ElementNode> *distribution_factors_field = NULL;
00487     bool surface_df_defined = false; // Has the surface df field been defined yet?
00488 
00489 
00490     int block_count = sset->block_count();
00491     for (int i=0; i < block_count; i++) {
00492       Ioss::SideBlock *side_block = sset->get_block(i);
00493       if (stk::io::include_entity(side_block)) {
00494   stk::mesh::Part * const side_block_part = meta.get_part(side_block->name());
00495   assert(side_block_part != NULL);
00496   meta.declare_part_subset(*fs_part, *side_block_part);
00497 
00498         const stk::mesh::EntityRank part_rank = side_block_part->primary_entity_rank();
00499 
00500   if (side_block->field_exists("distribution_factors")) {
00501     if (!surface_df_defined) {
00502       std::string field_name = sset->name() + "_distribution_factors";
00503       distribution_factors_field =
00504         &meta.declare_field<stk::mesh::Field<double, stk::mesh::ElementNode> >(field_name);
00505       stk::io::set_distribution_factor_field(*fs_part, *distribution_factors_field);
00506       surface_df_defined = true;
00507     }
00508     stk::io::set_distribution_factor_field(*side_block_part, *distribution_factors_field);
00509     int side_node_count = side_block->topology()->number_nodes();
00510     stk::mesh::put_field(*distribution_factors_field,
00511                                part_rank,
00512                                *side_block_part, side_node_count);
00513   }
00514 
00519   stk::io::define_io_fields(side_block, Ioss::Field::TRANSIENT,
00520           *side_block_part,
00521                                   part_rank);
00522       }
00523     }
00524   }
00525 
00526   // ========================================================================
00527   void process_sidesets(Ioss::Region &region, stk::mesh::MetaData &meta)
00528   {
00529     stk::mesh::fem::FEMMetaData &fem = stk::mesh::fem::FEMMetaData::get(meta);
00530     const stk::mesh::EntityRank side_rank = fem.side_rank();
00531 
00532     const Ioss::SideSetContainer& side_sets = region.get_sidesets();
00533     stk::io::default_part_processing(side_sets, meta, side_rank);
00534 
00535     for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
00536   it != side_sets.end(); ++it) {
00537       Ioss::SideSet *entity = *it;
00538 
00539       if (stk::io::include_entity(entity)) {
00540   process_surface_entity(entity, meta, side_rank);
00541       }
00542     }
00543   }
00544 
00545   // ========================================================================
00546   // Bulk Data
00547   // ========================================================================
00548   void process_nodeblocks(Ioss::Region &region, stk::mesh::BulkData &bulk)
00549   {
00550     // This must be called after the "process_element_blocks" call
00551     // since there may be nodes that exist in the database that are
00552     // not part of the analysis mesh due to subsetting of the element
00553     // blocks.
00554 
00555     const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
00556     assert(node_blocks.size() == 1);
00557 
00558     Ioss::NodeBlock *nb = node_blocks[0];
00559 
00560     std::vector<stk::mesh::Entity*> nodes;
00561     stk::io::get_entity_list(nb, stk::mesh::fem::FEMMetaData::NODE_RANK, bulk, nodes);
00562 
00567     const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk);
00568     stk::mesh::Field<double,stk::mesh::Cartesian> *coord_field =
00569       meta.get_field<stk::mesh::Field<double,stk::mesh::Cartesian> >("coordinates");
00570 
00571     stk::io::field_data_from_ioss(coord_field, nodes, nb, "mesh_model_coordinates");
00572   }
00573 
00574   // ========================================================================
00575   void process_elementblocks(Ioss::Region &region, stk::mesh::BulkData &bulk)
00576   {
00577     const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
00578 
00579     for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
00580   it != elem_blocks.end(); ++it) {
00581       Ioss::ElementBlock *entity = *it;
00582 
00583       if (stk::io::include_entity(entity)) {
00584   const std::string &name = entity->name();
00585         const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk);
00586         const stk::mesh::fem::FEMMetaData &fem_meta_data = stk::mesh::fem::FEMMetaData::get(meta);
00587   stk::mesh::Part* const part = meta.get_part(name);
00588   assert(part != NULL);
00589 
00590   const CellTopologyData* cell_topo = fem_meta_data.get_cell_topology(*part).getCellTopologyData();
00591   if (cell_topo == NULL) {
00592           std::ostringstream msg ;
00593     msg << " INTERNAL_ERROR: Part " << part->name() << " returned NULL from get_cell_topology()";
00594     throw std::runtime_error( msg.str() );
00595   }
00596 
00597   std::vector<int> elem_ids ;
00598   std::vector<int> connectivity ;
00599   std::vector<stk::mesh::EntityId> connectivity2 ;
00600 
00601   entity->get_field_data("ids", elem_ids);
00602   entity->get_field_data("connectivity", connectivity);
00603         connectivity2.reserve(connectivity.size());
00604         std::copy(connectivity.begin(), connectivity.end(), std::back_inserter(connectivity2));
00605 
00606   size_t element_count = elem_ids.size();
00607   int nodes_per_elem = cell_topo->node_count ;
00608 
00609   std::vector<stk::mesh::Entity*> elements(element_count);
00610   for(size_t i=0; i<element_count; ++i) {
00611     stk::mesh::EntityId *conn = &connectivity2[i*nodes_per_elem];
00612     elements[i] = &stk::mesh::fem::declare_element(bulk, *part, elem_ids[i], conn);
00613   }
00614 
00615   // For this example, we are just taking all attribute fields
00616   // found on the io database and populating fields on the
00617   // corresponding mesh part.  In practice, would probably be
00618   // selective about which attributes to use...
00619   Ioss::NameList names;
00620   entity->field_describe(Ioss::Field::ATTRIBUTE, &names);
00621   for (Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00622     if (*I == "attribute" && names.size() > 1)
00623       continue;
00624     stk::mesh::FieldBase *field = meta.get_field<stk::mesh::FieldBase>(*I);
00625     stk::io::field_data_from_ioss(field, elements, entity, *I);
00626 
00627   }
00628       }
00629     }
00630   }
00631 
00632   // ========================================================================
00633   void process_nodesets(Ioss::Region &region, stk::mesh::BulkData &bulk)
00634   {
00635     // Should only process nodes that have already been defined via the element
00636     // blocks connectivity lists.
00637     const Ioss::NodeSetContainer& node_sets = region.get_nodesets();
00638 
00639     for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
00640   it != node_sets.end(); ++it) {
00641       Ioss::NodeSet *entity = *it;
00642 
00643       if (stk::io::include_entity(entity)) {
00644   const std::string & name = entity->name();
00645   const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk);
00646   stk::mesh::Part* const part = meta.get_part(name);
00647   assert(part != NULL);
00648   stk::mesh::PartVector add_parts( 1 , part );
00649 
00650   std::vector<int> node_ids ;
00651   int node_count = entity->get_field_data("ids", node_ids);
00652 
00653   std::vector<stk::mesh::Entity*> nodes(node_count);
00654   for(int i=0; i<node_count; ++i) {
00655     nodes[i] = bulk.get_entity( stk::mesh::fem::FEMMetaData::NODE_RANK, node_ids[i] );
00656     if (nodes[i] != NULL)
00657       bulk.declare_entity(stk::mesh::fem::FEMMetaData::NODE_RANK, node_ids[i], add_parts );
00658   }
00659 
00664   stk::mesh::Field<double> *df_field =
00665     meta.get_field<stk::mesh::Field<double> >("distribution_factors");
00666 
00667   if (df_field != NULL) {
00668     stk::io::field_data_from_ioss(df_field, nodes, entity, "distribution_factors");
00669   }
00670       }
00671     }
00672   }
00673 
00674   // ========================================================================
00675   void process_surface_entity(const Ioss::SideSet* sset ,
00676             stk::mesh::BulkData & bulk)
00677   {
00678     assert(sset->type() == Ioss::SIDESET);
00679 
00680     const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk);
00681     const stk::mesh::fem::FEMMetaData &fem_meta_data = stk::mesh::fem::FEMMetaData::get(meta);
00682     const stk::mesh::EntityRank element_rank = fem_meta_data.element_rank();
00683 
00684     int block_count = sset->block_count();
00685     for (int i=0; i < block_count; i++) {
00686       Ioss::SideBlock *block = sset->get_block(i);
00687       if (stk::io::include_entity(block)) {
00688   std::vector<int> side_ids ;
00689   std::vector<int> elem_side ;
00690 
00691   stk::mesh::Part * const side_block_part = meta.get_part(block->name());
00692   stk::mesh::EntityRank side_rank = side_block_part->primary_entity_rank();
00693 
00694   block->get_field_data("ids", side_ids);
00695   block->get_field_data("element_side", elem_side);
00696 
00697   assert(side_ids.size() * 2 == elem_side.size());
00698   stk::mesh::PartVector add_parts( 1 , side_block_part );
00699 
00700   size_t side_count = side_ids.size();
00701   std::vector<stk::mesh::Entity*> sides(side_count);
00702   for(size_t is=0; is<side_count; ++is) {
00703 
00704     stk::mesh::Entity* const elem = bulk.get_entity(element_rank, elem_side[is*2]);
00705 
00706     // If NULL, then the element was probably assigned to an
00707     // element block that appears in the database, but was
00708     // subsetted out of the analysis mesh. Only process if
00709     // non-null.
00710     if (elem != NULL) {
00711       // Ioss uses 1-based side ordinal, stk::mesh uses 0-based.
00712       // Hence the '-1' in the following line.
00713       int side_ordinal = elem_side[is*2+1] - 1 ;
00714 
00715       stk::mesh::Entity *side = NULL;
00716       if (side_rank == 2) {
00717         side = &stk::mesh::fem::declare_element_side(bulk, side_ids[is], *elem, side_ordinal);
00718       } else {
00719         side = &stk::mesh::fem::declare_element_edge(bulk, side_ids[is], *elem, side_ordinal);
00720       }
00721       bulk.change_entity_parts( *side, add_parts );
00722       sides[is] = side;
00723     } else {
00724       sides[is] = NULL;
00725     }
00726   }
00727 
00728   const stk::mesh::Field<double, stk::mesh::ElementNode> *df_field =
00729     stk::io::get_distribution_factor_field(*side_block_part);
00730   if (df_field != NULL) {
00731     stk::io::field_data_from_ioss(df_field, sides, block, "distribution_factors");
00732   }
00733       }
00734     }
00735   }
00736 
00737   // ========================================================================
00738   void process_sidesets(Ioss::Region &region, stk::mesh::BulkData &bulk)
00739   {
00740     const Ioss::SideSetContainer& side_sets = region.get_sidesets();
00741 
00742     for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
00743   it != side_sets.end(); ++it) {
00744       Ioss::SideSet *entity = *it;
00745 
00746       if (stk::io::include_entity(entity)) {
00747   process_surface_entity(entity, bulk);
00748       }
00749     }
00750   }
00751 
00752   // ========================================================================
00753   // ========================================================================
00754   void get_field_data(stk::mesh::BulkData &bulk, stk::mesh::Part &part,
00755           stk::mesh::EntityRank part_type,
00756           Ioss::GroupingEntity *io_entity,
00757           Ioss::Field::RoleType filter_role)
00758   {
00759     std::vector<stk::mesh::Entity*> entities;
00760     stk::io::get_entity_list(io_entity, part_type, bulk, entities);
00761 
00762     stk::mesh::MetaData& meta = stk::mesh::MetaData::get(part);
00763     stk::mesh::Part &universal = meta.universal_part();
00764     const std::vector<stk::mesh::FieldBase*> &fields = meta.get_fields();
00765 
00766     std::vector<stk::mesh::FieldBase *>::const_iterator I = fields.begin();
00767     while (I != fields.end()) {
00768       const stk::mesh::FieldBase *f = *I; ++I;
00769       if (stk::io::is_valid_part_field(f, part_type, part, universal, filter_role)) {
00770   stk::io::field_data_from_ioss(f, entities, io_entity, f->name());
00771       }
00772     }
00773   }
00774 
00775   void process_input_request(Ioss::Region &region,
00776            stk::mesh::BulkData &bulk,
00777            int step)
00778   {
00779     region.begin_state(step);
00780 
00781     // Special processing for nodeblock (all nodes in model)...
00782     const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk);
00783 
00784     // ??? Get field data from nodeblock...
00785     get_field_data(bulk, meta.universal_part(), stk::mesh::fem::FEMMetaData::NODE_RANK,
00786        region.get_node_blocks()[0], Ioss::Field::TRANSIENT);
00787 
00788     const stk::mesh::PartVector & all_parts = meta.get_parts();
00789     for ( stk::mesh::PartVector::const_iterator
00790       ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
00791 
00792       stk::mesh::Part * const part = *ip;
00793 
00794       const stk::mesh::EntityRank part_rank = part->primary_entity_rank();
00795 
00796       // Check whether this part should be output to results database.
00797       if (stk::io::is_part_io_part(*part)) {
00798   // Get Ioss::GroupingEntity corresponding to this part...
00799   Ioss::GroupingEntity *entity = region.get_entity(part->name());
00800   if (entity != NULL) {
00801     if (entity->type() == Ioss::SIDESET) {
00802       Ioss::SideSet *sset = dynamic_cast<Ioss::SideSet*>(entity);
00803       assert(sset != NULL);
00804       int block_count = sset->block_count();
00805       for (int i=0; i < block_count; i++) {
00806         Ioss::SideBlock *side_block = sset->get_block(i);
00808     get_field_data(bulk, *part,
00809                                part_rank,
00810              side_block, Ioss::Field::TRANSIENT);
00811       }
00812     } else {
00813       get_field_data(bulk, *part,
00814                            part_rank,
00815          entity, Ioss::Field::TRANSIENT);
00816     }
00817   } else {
00820   }
00821       }
00822     }
00823 
00824     region.end_state(step);
00825   }
00826 
00827   void put_field_data(stk::mesh::BulkData &bulk, stk::mesh::Part &part,
00828           stk::mesh::EntityRank part_type,
00829           Ioss::GroupingEntity *io_entity,
00830           Ioss::Field::RoleType filter_role)
00831   {
00832     std::vector<stk::mesh::Entity*> entities;
00833     stk::io::get_entity_list(io_entity, part_type, bulk, entities);
00834 
00835     stk::mesh::MetaData& meta = stk::mesh::MetaData::get(part);
00836     stk::mesh::Part &universal = meta.universal_part();
00837     const std::vector<stk::mesh::FieldBase*> &fields = meta.get_fields();
00838 
00839     std::vector<stk::mesh::FieldBase *>::const_iterator I = fields.begin();
00840     while (I != fields.end()) {
00841       const stk::mesh::FieldBase *f = *I; ++I;
00842       if (stk::io::is_valid_part_field(f, part_type, part, universal, filter_role)) {
00843   stk::io::field_data_to_ioss(f, entities, io_entity, f->name(), filter_role);
00844       }
00845     }
00846   }
00847 
00848   void process_output_request(Ioss::Region &region,
00849             stk::mesh::BulkData &bulk,
00850             int step)
00851   {
00852     region.begin_state(step);
00853     // Special processing for nodeblock (all nodes in model)...
00854     const stk::mesh::MetaData& meta = stk::mesh::MetaData::get(bulk);
00855 
00856     put_field_data(bulk, meta.universal_part(), stk::mesh::fem::FEMMetaData::NODE_RANK,
00857        region.get_node_blocks()[0], Ioss::Field::TRANSIENT);
00858 
00859     const stk::mesh::PartVector & all_parts = meta.get_parts();
00860     for ( stk::mesh::PartVector::const_iterator
00861       ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
00862 
00863       stk::mesh::Part * const part = *ip;
00864 
00865       const stk::mesh::EntityRank part_rank = part->primary_entity_rank();
00866 
00867       // Check whether this part should be output to results database.
00868       if (stk::io::is_part_io_part(*part)) {
00869 
00870   // Get Ioss::GroupingEntity corresponding to this part...
00871   Ioss::GroupingEntity *entity = region.get_entity(part->name());
00872   if (entity != NULL) {
00873 
00874     if (entity->type() == Ioss::SIDESET) {
00875       Ioss::SideSet *sset = dynamic_cast<Ioss::SideSet*>(entity);
00876       assert(sset != NULL);
00877       int block_count = sset->block_count();
00878 
00879       for (int i=0; i < block_count; i++) {
00880         Ioss::SideBlock *side_block = sset->get_block(i);
00882     put_field_data(bulk, *part, part_rank,
00883              side_block, Ioss::Field::TRANSIENT);
00884       }
00885     } else {
00886       put_field_data(bulk, *part, part_rank,
00887          entity, Ioss::Field::TRANSIENT);
00888     }
00889   } else {
00892   }
00893       }
00894     }
00895     region.end_state(step);
00896   }
00897 }
00898 
00899 // ========================================================================
00900 #include <boost/program_options.hpp>
00901 
00902 #include <stk_util/parallel/BroadcastArg.hpp>
00903 #include <stk_util/environment/ProgramOptions.hpp>
00904 
00905   namespace bopt = boost::program_options;
00906   int main(int argc, char** argv)
00907   {
00908     //----------------------------------
00909     // Broadcast argc and argv to all processors.
00910 
00911     stk::ParallelMachine comm = stk::parallel_machine_init(&argc, &argv);
00912 
00913     stk::BroadcastArg b_arg(comm, argc, argv);
00914 
00915     //----------------------------------
00916     // Process the broadcast command line arguments
00917 
00918     bopt::options_description desc("options");
00919 
00920     desc.add_options()
00921       ("help,h",        "produce help message")
00922       ("mesh",         bopt::value<std::string>(), "mesh file" )
00923       ("directory,d",  bopt::value<std::string>(), "working directory" )
00924       ("output-log,o", bopt::value<std::string>(), "output log path" )
00925       ("runtest,r",    bopt::value<std::string>(), "runtest pid file" );
00926 
00927     stk::get_options_description().add(desc);
00928 
00929     bopt::variables_map &vm = stk::get_variables_map();
00930     try {
00931       bopt::store(bopt::parse_command_line(b_arg.m_argc, b_arg.m_argv, desc), vm);
00932       bopt::notify(vm);
00933     }
00934     catch (std::exception & /* x */) {
00935       std::exit(1);
00936     }
00937 
00938     if (vm.count("help")) {
00939       std::cout << desc << "\n";
00940       std::exit(EXIT_SUCCESS);
00941     }
00942 
00943     //----------------------------------
00944 
00945     if ( vm.count("mesh") ) {
00946       std::string in_filename = boost::any_cast<std::string>(vm["mesh"].value());
00947       std::string out_filename = in_filename + ".out";
00948       stk_example_io::io_example(comm, in_filename, out_filename );
00949     } else {
00950       std::cout << "OPTION ERROR: The '--mesh <filename>' option is required!\n";
00951       std::exit(EXIT_FAILURE);
00952     }
00953     stk::parallel_machine_finalize();
00954 
00955     return 0;
00956   }
00957 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends