Sierra Toolkit Version of the Day
MeshReadWriteUtils.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 <stk_io/MeshReadWriteUtils.hpp>
00010 
00011 #include <stk_mesh/base/BulkData.hpp>
00012 #include <stk_mesh/base/GetEntities.hpp>
00013 #include <stk_mesh/fem/FEMMetaData.hpp>
00014 #include <stk_mesh/fem/FEMHelpers.hpp>
00015 
00016 #include <stk_mesh/base/Field.hpp>
00017 #include <stk_mesh/base/FieldData.hpp>
00018 #include <stk_mesh/base/FieldParallel.hpp>
00019 
00020 #include <Shards_BasicTopologies.hpp>
00021 
00022 #include <stk_io/IossBridge.hpp>
00023 
00024 #include <Ioss_SubSystem.h>
00025 
00026 #include <stk_util/util/tokenize.hpp>
00027 #include <iostream>
00028 #include <sstream>
00029 #include <cmath>
00030 
00031 #include <limits>
00032 #include <assert.h>
00033 
00034 namespace {
00035 void process_surface_entity(Ioss::SideSet *sset, stk::mesh::fem::FEMMetaData &fem_meta)
00036 {
00037   assert(sset->type() == Ioss::SIDESET);
00038   const Ioss::SideBlockContainer& blocks = sset->get_side_blocks();
00039   stk::io::default_part_processing(blocks, fem_meta);
00040 
00041   stk::mesh::Part* const ss_part = fem_meta.get_part(sset->name());
00042   assert(ss_part != NULL);
00043 
00044   stk::mesh::Field<double, stk::mesh::ElementNode> *distribution_factors_field = NULL;
00045   bool surface_df_defined = false; // Has the surface df field been defined yet?
00046 
00047   size_t block_count = sset->block_count();
00048   for (size_t i=0; i < block_count; i++) {
00049     Ioss::SideBlock *sb = sset->get_block(i);
00050     if (stk::io::include_entity(sb)) {
00051       stk::mesh::Part * const sb_part = fem_meta.get_part(sb->name());
00052       assert(sb_part != NULL);
00053       fem_meta.declare_part_subset(*ss_part, *sb_part);
00054 
00055       if (sb->field_exists("distribution_factors")) {
00056         if (!surface_df_defined) {
00057           std::string field_name = sset->name() + "_df";
00058           distribution_factors_field =
00059             &fem_meta.declare_field<stk::mesh::Field<double, stk::mesh::ElementNode> >(field_name);
00060           stk::io::set_field_role(*distribution_factors_field, Ioss::Field::MESH);
00061           stk::io::set_distribution_factor_field(*ss_part, *distribution_factors_field);
00062           surface_df_defined = true;
00063         }
00064         stk::io::set_distribution_factor_field(*sb_part, *distribution_factors_field);
00065         int side_node_count = sb->topology()->number_nodes();
00066         stk::mesh::put_field(*distribution_factors_field,
00067                              stk::io::part_primary_entity_rank(*sb_part),
00068                              *sb_part, side_node_count);
00069       }
00070     }
00071   }
00072 }
00073 
00074   size_t get_entities(stk::mesh::Part &part,
00075           const stk::mesh::BulkData &bulk,
00076           std::vector<stk::mesh::Entity*> &entities,
00077           const stk::mesh::Selector *anded_selector)
00078   {
00079     stk::mesh::MetaData & meta = stk::mesh::MetaData::get(part);
00080     stk::mesh::EntityRank type = stk::io::part_primary_entity_rank(part);
00081     
00082     stk::mesh::Selector own = meta.locally_owned_part();
00083     stk::mesh::Selector selector = part & own;
00084     if (anded_selector) selector &= *anded_selector;
00085     
00086     get_selected_entities(selector, bulk.buckets(type), entities);
00087     return entities.size();
00088   }
00089 }
00090 
00091 // ========================================================================
00092 template <typename INT>
00093 void process_surface_entity(const Ioss::SideSet* sset, stk::mesh::BulkData & bulk, INT /*dummy*/)
00094 {
00095   assert(sset->type() == Ioss::SIDESET);
00096 
00097   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00098 
00099   size_t block_count = sset->block_count();
00100   for (size_t i=0; i < block_count; i++) {
00101     Ioss::SideBlock *block = sset->get_block(i);
00102     if (stk::io::include_entity(block)) {
00103       std::vector<INT> side_ids ;
00104       std::vector<INT> elem_side ;
00105 
00106       stk::mesh::Part * const sb_part = fem_meta.get_part(block->name());
00107       stk::mesh::EntityRank elem_rank = fem_meta.element_rank();
00108 
00109       block->get_field_data("ids", side_ids);
00110       block->get_field_data("element_side", elem_side);
00111 
00112       assert(side_ids.size() * 2 == elem_side.size());
00113       stk::mesh::PartVector add_parts( 1 , sb_part );
00114 
00115       size_t side_count = side_ids.size();
00116       std::vector<stk::mesh::Entity*> sides(side_count);
00117       for(size_t is=0; is<side_count; ++is) {
00118         stk::mesh::Entity* const elem = bulk.get_entity(elem_rank, elem_side[is*2]);
00119 
00120         // If NULL, then the element was probably assigned to an
00121         // element block that appears in the database, but was
00122         // subsetted out of the analysis mesh. Only process if
00123         // non-null.
00124         if (elem != NULL) {
00125           // Ioss uses 1-based side ordinal, stk::mesh uses 0-based.
00126           int side_ordinal = elem_side[is*2+1] - 1;
00127 
00128           stk::mesh::Entity* side_ptr = NULL;
00129           side_ptr = &stk::mesh::fem::declare_element_side(bulk, side_ids[is], *elem, side_ordinal);
00130           stk::mesh::Entity& side = *side_ptr;
00131 
00132           bulk.change_entity_parts( side, add_parts );
00133           sides[is] = &side;
00134         } else {
00135           sides[is] = NULL;
00136         }
00137       }
00138 
00139       const stk::mesh::Field<double, stk::mesh::ElementNode> *df_field =
00140         stk::io::get_distribution_factor_field(*sb_part);
00141       if (df_field != NULL) {
00142         stk::io::field_data_from_ioss(df_field, sides, block, "distribution_factors");
00143       }
00144 
00145       // Add all attributes as fields.
00146       // If the only attribute is 'attribute', then add it; otherwise the other attributes are the
00147       // named components of the 'attribute' field, so add them instead.
00148       Ioss::NameList names;
00149       block->field_describe(Ioss::Field::ATTRIBUTE, &names);
00150       for(Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00151         if(*I == "attribute" && names.size() > 1)
00152           continue;
00153         stk::mesh::FieldBase *field = fem_meta.get_field<stk::mesh::FieldBase> (*I);
00154         if (field)
00155           stk::io::field_data_from_ioss(field, sides, block, *I);
00156       }
00157     }
00158   }
00159 }
00160 
00161 void process_surface_entity(const Ioss::SideSet* sset, stk::mesh::BulkData & bulk)
00162 {
00163   if (stk::io::db_api_int_size(sset) == 4)
00164     process_surface_entity(sset, bulk, (int)0);
00165   else
00166     process_surface_entity(sset, bulk, (int64_t)0);
00167 }
00168 
00169 void process_nodeblocks(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00170 {
00171   const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
00172   assert(node_blocks.size() == 1);
00173 
00174   Ioss::NodeBlock *nb = node_blocks[0];
00175 
00176   assert(nb->field_exists("mesh_model_coordinates"));
00177   Ioss::Field coordinates = nb->get_field("mesh_model_coordinates");
00178   int spatial_dim = coordinates.transformed_storage()->component_count();
00179 
00180   stk::mesh::Field<double,stk::mesh::Cartesian> & coord_field =
00181     fem_meta.declare_field<stk::mesh::Field<double,stk::mesh::Cartesian> >("coordinates");
00182 
00183   stk::mesh::put_field( coord_field, fem_meta.node_rank(), fem_meta.universal_part(), spatial_dim);
00184   stk::io::define_io_fields(nb, Ioss::Field::ATTRIBUTE, fem_meta.universal_part(), 0);
00185 }
00186 
00187 void process_nodeblocks(Ioss::Region &region, stk::mesh::BulkData &bulk)
00188 {
00189   // This must be called after the "process_element_blocks" call
00190   // since there may be nodes that exist in the database that are
00191   // not part of the analysis mesh due to subsetting of the element
00192   // blocks.
00193 
00194   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00195 
00196   const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
00197   assert(node_blocks.size() == 1);
00198 
00199   Ioss::NodeBlock *nb = node_blocks[0];
00200 
00201   std::vector<stk::mesh::Entity*> nodes;
00202   stk::io::get_entity_list(nb, fem_meta.node_rank(), bulk, nodes);
00203 
00204   stk::mesh::Field<double,stk::mesh::Cartesian> *coord_field =
00205     fem_meta.get_field<stk::mesh::Field<double,stk::mesh::Cartesian> >("coordinates");
00206 
00207   stk::io::field_data_from_ioss(coord_field, nodes, nb, "mesh_model_coordinates");
00208 
00209   // Add all attributes as fields.
00210   // If the only attribute is 'attribute', then add it; otherwise the other attributes are the
00211   // named components of the 'attribute' field, so add them instead.
00212   Ioss::NameList names;
00213   nb->field_describe(Ioss::Field::ATTRIBUTE, &names);
00214   for(Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00215     if(*I == "attribute" && names.size() > 1)
00216       continue;
00217     stk::mesh::FieldBase *field = fem_meta.get_field<stk::mesh::FieldBase> (*I);
00218     if (field)
00219       stk::io::field_data_from_ioss(field, nodes, nb, *I);
00220   }
00221 }
00222 
00223 // ========================================================================
00224 void process_elementblocks(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00225 {
00226   const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
00227   stk::io::default_part_processing(elem_blocks, fem_meta);
00228 }
00229 
00230 template <typename INT>
00231 void process_elementblocks(Ioss::Region &region, stk::mesh::BulkData &bulk, INT /*dummy*/)
00232 {
00233   const stk::mesh::fem::FEMMetaData& fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00234 
00235   const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
00236   for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
00237       it != elem_blocks.end(); ++it) {
00238     Ioss::ElementBlock *entity = *it;
00239 
00240     if (stk::io::include_entity(entity)) {
00241       const std::string &name = entity->name();
00242       stk::mesh::Part* const part = fem_meta.get_part(name);
00243       assert(part != NULL);
00244 
00245       const CellTopologyData* cell_topo = stk::io::get_cell_topology(*part);
00246       if (cell_topo == NULL) {
00247         std::ostringstream msg ;
00248         msg << " INTERNAL_ERROR: Part " << part->name() << " returned NULL from get_cell_topology()";
00249         throw std::runtime_error( msg.str() );
00250       }
00251 
00252       std::vector<INT> elem_ids ;
00253       std::vector<INT> connectivity ;
00254 
00255       entity->get_field_data("ids", elem_ids);
00256       entity->get_field_data("connectivity", connectivity);
00257 
00258       size_t element_count = elem_ids.size();
00259       int nodes_per_elem = cell_topo->node_count ;
00260 
00261       std::vector<stk::mesh::EntityId> id_vec(nodes_per_elem);
00262       std::vector<stk::mesh::Entity*> elements(element_count);
00263 
00264       for(size_t i=0; i<element_count; ++i) {
00265         INT *conn = &connectivity[i*nodes_per_elem];
00266         std::copy(&conn[0], &conn[0+nodes_per_elem], id_vec.begin());
00267         elements[i] = &stk::mesh::fem::declare_element(bulk, *part, elem_ids[i], &id_vec[0]);
00268       }
00269 
00270       // Add all element attributes as fields.
00271       // If the only attribute is 'attribute', then add it; otherwise the other attributes are the
00272       // named components of the 'attribute' field, so add them instead.
00273       Ioss::NameList names;
00274       entity->field_describe(Ioss::Field::ATTRIBUTE, &names);
00275       for(Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00276         if(*I == "attribute" && names.size() > 1)
00277           continue;
00278         stk::mesh::FieldBase *field = fem_meta.get_field<stk::mesh::FieldBase> (*I);
00279         if (field)
00280           stk::io::field_data_from_ioss(field, elements, entity, *I);
00281       }
00282     }
00283   }
00284 }
00285 
00286 // ========================================================================
00287 // ========================================================================
00288 void process_nodesets(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00289 {
00290   const Ioss::NodeSetContainer& node_sets = region.get_nodesets();
00291   stk::io::default_part_processing(node_sets, fem_meta);
00292 
00293   stk::mesh::Field<double> & distribution_factors_field =
00294     fem_meta.declare_field<stk::mesh::Field<double> >("distribution_factors");
00295   stk::io::set_field_role(distribution_factors_field, Ioss::Field::MESH);
00296 
00302   for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
00303       it != node_sets.end(); ++it) {
00304     Ioss::NodeSet *entity = *it;
00305 
00306     if (stk::io::include_entity(entity)) {
00307       stk::mesh::Part* const part = fem_meta.get_part(entity->name());
00308       assert(part != NULL);
00309       assert(entity->field_exists("distribution_factors"));
00310 
00311       stk::mesh::put_field(distribution_factors_field, fem_meta.node_rank(), *part);
00312     }
00313   }
00314 }
00315 
00316 // ========================================================================
00317 // ========================================================================
00318 void process_sidesets(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00319 {
00320   const Ioss::SideSetContainer& side_sets = region.get_sidesets();
00321   stk::io::default_part_processing(side_sets, fem_meta);
00322 
00323   for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
00324       it != side_sets.end(); ++it) {
00325     Ioss::SideSet *entity = *it;
00326 
00327     if (stk::io::include_entity(entity)) {
00328       process_surface_entity(entity, fem_meta);
00329     }
00330   }
00331 }
00332 
00333 // ========================================================================
00334 template <typename INT>
00335 void process_nodesets(Ioss::Region &region, stk::mesh::BulkData &bulk, INT /*dummy*/)
00336 {
00337   // Should only process nodes that have already been defined via the element
00338   // blocks connectivity lists.
00339   const Ioss::NodeSetContainer& node_sets = region.get_nodesets();
00340   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00341 
00342   for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
00343       it != node_sets.end(); ++it) {
00344     Ioss::NodeSet *entity = *it;
00345 
00346     if (stk::io::include_entity(entity)) {
00347       const std::string & name = entity->name();
00348       stk::mesh::Part* const part = fem_meta.get_part(name);
00349       assert(part != NULL);
00350       stk::mesh::PartVector add_parts( 1 , part );
00351 
00352       std::vector<INT> node_ids ;
00353       size_t node_count = entity->get_field_data("ids", node_ids);
00354 
00355       std::vector<stk::mesh::Entity*> nodes(node_count);
00356       stk::mesh::EntityRank n_rank = fem_meta.node_rank();
00357       for(size_t i=0; i<node_count; ++i) {
00358         nodes[i] = bulk.get_entity(n_rank, node_ids[i] );
00359         if (nodes[i] != NULL)
00360           bulk.declare_entity(n_rank, node_ids[i], add_parts );
00361       }
00362 
00363       stk::mesh::Field<double> *df_field =
00364         fem_meta.get_field<stk::mesh::Field<double> >("distribution_factors");
00365 
00366       if (df_field != NULL) {
00367         stk::io::field_data_from_ioss(df_field, nodes, entity, "distribution_factors");
00368       }
00369 
00370       // Add all attributes as fields.
00371       // If the only attribute is 'attribute', then add it; otherwise the other attributes are the
00372       // named components of the 'attribute' field, so add them instead.
00373       Ioss::NameList names;
00374       entity->field_describe(Ioss::Field::ATTRIBUTE, &names);
00375       for(Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00376         if(*I == "attribute" && names.size() > 1)
00377           continue;
00378         stk::mesh::FieldBase *field = fem_meta.get_field<stk::mesh::FieldBase> (*I);
00379         if (field)
00380           stk::io::field_data_from_ioss(field, nodes, entity, *I);
00381       }
00382     }
00383   }
00384 }
00385 
00386 // ========================================================================
00387 void process_sidesets(Ioss::Region &region, stk::mesh::BulkData &bulk)
00388 {
00389   const Ioss::SideSetContainer& side_sets = region.get_sidesets();
00390 
00391   for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
00392       it != side_sets.end(); ++it) {
00393     Ioss::SideSet *entity = *it;
00394 
00395     if (stk::io::include_entity(entity)) {
00396       process_surface_entity(entity, bulk);
00397     }
00398   }
00399 }
00400 
00401 // ========================================================================
00402 void put_field_data(stk::mesh::BulkData &bulk, stk::mesh::Part &part,
00403                     stk::mesh::EntityRank part_type,
00404                     Ioss::GroupingEntity *io_entity,
00405                     Ioss::Field::RoleType filter_role,
00406         const stk::mesh::Selector *anded_selector=NULL)
00407 {
00408   std::vector<stk::mesh::Entity*> entities;
00409   if (io_entity->type() == Ioss::SIDEBLOCK) {
00410     // Temporary Kluge to handle sideblocks which contain internally generated sides
00411     // where the "ids" field on the io_entity doesn't work to get the correct side...
00412     // NOTE: Could use this method for all entity types, but then need to correctly 
00413     // specify whether shared entities are included/excluded (See IossBridge version).
00414     size_t num_sides = get_entities(part, bulk, entities, anded_selector);
00415     if (num_sides != (size_t)io_entity->get_property("entity_count").get_int()) {
00416       std::ostringstream msg ;
00417       msg << " INTERNAL_ERROR: Number of sides on part " << part.name() << " (" << num_sides
00418     << ") does not match number of sides in the associated Ioss SideBlock named "
00419     << io_entity->name() << " (" << io_entity->get_property("entity_count").get_int()
00420     << ").";
00421         throw std::runtime_error( msg.str() );
00422     }
00423   } else {
00424     stk::io::get_entity_list(io_entity, part_type, bulk, entities);
00425   }
00426 
00427   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00428   const std::vector<stk::mesh::FieldBase*> &fields = fem_meta.get_fields();
00429 
00430   std::vector<stk::mesh::FieldBase *>::const_iterator I = fields.begin();
00431   while (I != fields.end()) {
00432     const stk::mesh::FieldBase *f = *I; ++I;
00433     stk::io::field_data_to_ioss(f, entities, io_entity, f->name(), filter_role);
00434   }
00435 }
00436 
00437 void internal_process_output_request(stk::io::MeshData &mesh_data,
00438                                      stk::mesh::BulkData &bulk,
00439                                      int step,
00440                                      const std::set<const stk::mesh::Part*> &exclude)
00441 {
00442   Ioss::Region *region = mesh_data.m_output_region;
00443   region->begin_state(step);
00444   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00445 
00446   // Special processing for nodeblock (all nodes in model)...
00447   put_field_data(bulk, fem_meta.universal_part(), fem_meta.node_rank(),
00448                  region->get_node_blocks()[0], Ioss::Field::Field::TRANSIENT,
00449      mesh_data.m_anded_selector);
00450 
00451   // Now handle all non-nodeblock parts...
00452   const stk::mesh::PartVector & all_parts = fem_meta.get_parts();
00453   for ( stk::mesh::PartVector::const_iterator
00454           ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
00455 
00456     stk::mesh::Part * const part = *ip;
00457 
00458     // Check whether this part should be output to results database.
00459     if (stk::io::is_part_io_part(*part) && !exclude.count(part)) {
00460       // Get Ioss::GroupingEntity corresponding to this part...
00461       Ioss::GroupingEntity *entity = region->get_entity(part->name());
00462       if (entity != NULL && entity->type() != Ioss::SIDESET) {
00463         put_field_data(bulk, *part, stk::io::part_primary_entity_rank(*part),
00464                        entity, Ioss::Field::Field::TRANSIENT,
00465            mesh_data.m_anded_selector);
00466       }
00467     }
00468   }
00469   region->end_state(step);
00470 }
00471 
00472 namespace stk {
00473 namespace io {
00474 
00475 MeshData::~MeshData()
00476 {
00477   delete m_input_region;
00478   delete m_output_region;
00479 }
00480 
00481 void show_mesh_help()
00482 {
00483   std::cerr << "Options are:\n"
00484             << "\n"
00485             << "filename -- specify the name of the file from which to read the\n"
00486             << "            mesh file. If the --directory option is specified, it will be\n"
00487             << "            prepended to the filename unless the filename specifies an absolute path.\n"
00488             << "\n"
00489             << "gen:NxMxL -- internally generate a hex mesh of size N by M by L\n"
00490             << "             intervals. See 'Generated Options' below for more options.\n"
00491             << "\n"
00492             << "Generated Options:\n"
00493             << "shell:xXyYzZ\n"
00494             << "The argument specifies whether there is a shell block\n"
00495             << "at the location. 'x' is minX, 'X' is maxX, etc.\n"
00496             << "\n"
00497             << "help -- no argument, shows valid options\n"
00498             << "\n"
00499             << "show -- no argument, prints out a summary of the settings used to\n"
00500             << "generate the mesh. The output will look similar to:\n"
00501             << "    \"10x12x8|shell:xX|bbox:-10,-10,-10,10,10,10|show\"\n"
00502             << "\n"
00503             << "    Mesh Parameters:\n"
00504             << "\tIntervals: 10 by 12 by 8\n"
00505             << "\tX = 2       * (0..10) + -10     Range: -10 <= X <= 10\n"
00506             << "\tY = 1.66667 * (0..12) + -10     Range: -10 <= Y <= 10\n"
00507             << "\tZ = 2.5     * (0..8)  + -10     Range: -10 <= Z <= 10\n"
00508             << "\tNode Count (total)    = 1287\n"
00509             << "\tElement Count (total) = 1152\n"
00510             << "\tBlock Count           = 3\n"
00511             << "\n"
00512             << "shell:xXyYzZ \n"
00513             << "which specifies whether there is a shell block at that\n"
00514             << "location. 'x' is minimum x face, 'X' is maximum x face,\n"
00515             << "similarly for y and z.  Note that the argument string is a\n"
00516             << "single multicharacter string.  You can add multiple shell blocks\n"
00517             << "to a face, for example, shell:xxx would add three layered shell\n"
00518             << "blocks on the minimum x face.  An error is output if a non\n"
00519             << "xXyYzZ character is found, but execution continues.\n"
00520             << "\n"
00521             << "zdecomp:n0 n1,n2,...,n#proc-1\n"
00522             << "which are the number of intervals in the z direction for each\n"
00523             << "processor in a pallel run.  If this option is specified, then\n"
00524             << "the total number of intervals in the z direction is the sum of\n"
00525             << "the n0, n1, ... An interval count must be specified for each\n"
00526             << "processor.  If this option is not specified, then the number of\n"
00527             << "intervals on each processor in the z direction is numZ/numProc\n"
00528             << "with the extras added to the lower numbered processors.\n"
00529             << "\n"
00530             << "scale:xs,ys,zs\n"
00531             << "which are the scale factors in the x, y, and z directions. All\n"
00532             << "three must be specified if this option is present.\n"
00533             << "\n"
00534             << "- offset -- argument = xoff, yoff, zoff which are the offsets in the\n"
00535             << "x, y, and z directions.  All three must be specified if this option\n"
00536             << "is present.\n"
00537             << "\n"
00538             << "- bbox -- argument = xmin, ymin, zmin, xmax, ymax, zmax\n"
00539             << "which specify the lower left and upper right corners of\n"
00540             << "the bounding box for the generated mesh.  This will\n"
00541             << "calculate the scale and offset which will fit the mesh in\n"
00542             << "the specified box.  All calculations are based on the currently\n"
00543             << "active interval settings. If scale or offset or zdecomp\n"
00544             << "specified later in the option list, you may not get the\n"
00545             << "desired bounding box.\n"
00546             << "\n"
00547             << "- rotate -- argument = axis,angle,axis,angle,...\n"
00548             << "where axis is 'x', 'y', or 'z' and angle is the rotation angle in\n"
00549             << "degrees. Multiple rotations are cumulative. The composite rotation\n"
00550             << "matrix is applied at the time the coordinates are retrieved after\n"
00551             << "scaling and offset are applied.\n"
00552             << "\n"
00553             << "The unrotated coordinate of a node at grid location i,j,k is:\n"
00554             << "\n"
00555             << "\tx = x_scale * i + x_off,\n"
00556             << "\ty = z_scale * j + y_off,\n"
00557             << "\tz = z_scale * k + z_off,\n"
00558             << "\n"
00559             << "The extent of the unrotated mesh will be:\n"
00560             << "\n"
00561             << "\tx_off <= x <= x_scale * numX + x_off\n"
00562             << "\ty_off <= y <= y_scale * numY + y_off\n"
00563             << "\tz_off <= z <= z_scale * numZ + z_off\n"
00564             << "\n"
00565             << "If an unrecognized option is specified, an error message will be\n"
00566             << "output and execution will continue.\n"
00567             << "\n"
00568             << "An example of valid input is:\n"
00569             << "\n"
00570             << "\t\"10x20x40|scale:1,0.5,0.25|offset:-5,-5,-5|shell:xX\"\n"
00571             << "\n"
00572             << "\n"
00573             << "This would create a mesh with 10 intervals in x, 20 in y, 40 in z\n"
00574             << "The mesh would be centered on 0,0,0 with a range of 10 in each\n"
00575             << "direction. There would be a shell layer on the min and max\n"
00576             << "x faces.\n"
00577             << "\n"
00578             << "NOTE: All options are processed in the order they appear in\n"
00579             << "the parameters string (except rotate which is applied at the\n"
00580             << "time the coordinates are generated/retrieved)\n"
00581             << "\n";
00582 }
00583 
00584 void create_input_mesh(const std::string &mesh_type,
00585                        const std::string &mesh_filename,
00586                        stk::ParallelMachine comm,
00587                        stk::mesh::fem::FEMMetaData &fem_meta,
00588                        stk::io::MeshData &mesh_data)
00589 {
00590   Ioss::Region *in_region = mesh_data.m_input_region;
00591   if (in_region == NULL) {
00592   // If in_region is NULL, then open the file;
00593   // If in_region is non-NULL, then user has given us a valid Ioss::Region that
00594   // should be used.
00595   Ioss::DatabaseIO *dbi = Ioss::IOFactory::create(mesh_type, mesh_filename,
00596                                                     Ioss::READ_MODEL, comm,
00597                                                     mesh_data.m_property_manager);
00598   if (dbi == NULL || !dbi->ok()) {
00599     std::cerr  << "ERROR: Could not open database '" << mesh_filename
00600                  << "' of type '" << mesh_type << "'\n";
00601     Ioss::NameList db_types;
00602     Ioss::IOFactory::describe(&db_types);
00603     std::cerr << "\nSupported database types:\n\t";
00604     for (Ioss::NameList::const_iterator IF = db_types.begin(); IF != db_types.end(); ++IF) {
00605       std::cerr << *IF << "  ";
00606     }
00607     std::cerr << "\n\n";
00608   }
00609 
00610   // NOTE: 'in_region' owns 'dbi' pointer at this time...
00611   in_region = new Ioss::Region(dbi, "input_model");
00612   mesh_data.m_input_region = in_region;
00613   }
00614 
00615   size_t spatial_dimension = in_region->get_property("spatial_dimension").get_int();
00616   initialize_spatial_dimension(fem_meta, spatial_dimension, stk::mesh::fem::entity_rank_names(spatial_dimension));
00617 
00618   process_elementblocks(*in_region, fem_meta);
00619   process_nodeblocks(*in_region,    fem_meta);
00620   process_sidesets(*in_region,      fem_meta);
00621   process_nodesets(*in_region,      fem_meta);
00622 }
00623 
00624 
00625 void create_output_mesh(const std::string &filename,
00626                         stk::ParallelMachine comm,
00627                         stk::mesh::BulkData &bulk_data,
00628                         MeshData &mesh_data)
00629 {
00630   Ioss::Region *out_region = NULL;
00631 
00632   std::string out_filename = filename;
00633   if (filename.empty()) {
00634   out_filename = "default_output_mesh";
00635   } else {
00636   // These filenames may be coming from the generated options which
00637   // may have forms similar to: "2x2x1|size:.05|height:-0.1,1"
00638   // Strip the name at the first "+:|," character:
00639   std::vector<std::string> tokens;
00640   stk::util::tokenize(out_filename, "+|:,", tokens);
00641   out_filename = tokens[0];
00642   }
00643 
00644   Ioss::DatabaseIO *dbo = Ioss::IOFactory::create("exodusII", out_filename,
00645                                                   Ioss::WRITE_RESULTS,
00646                                                   comm, mesh_data.m_property_manager);
00647   if (dbo == NULL || !dbo->ok()) {
00648   std::cerr << "ERROR: Could not open results database '" << out_filename
00649               << "' of type 'exodusII'\n";
00650   std::exit(EXIT_FAILURE);
00651   }
00652 
00653   // NOTE: 'out_region' owns 'dbo' pointer at this time...
00654   out_region = new Ioss::Region(dbo, "results_output");
00655 
00656   stk::io::define_output_db(*out_region, bulk_data, mesh_data.m_input_region, mesh_data.m_anded_selector);
00657   stk::io::write_output_db(*out_region,  bulk_data, mesh_data.m_anded_selector);
00658   mesh_data.m_output_region = out_region;
00659 }
00660 
00661 // ========================================================================
00662 int process_output_request(MeshData &mesh_data,
00663                            stk::mesh::BulkData &bulk,
00664                            double time,
00665                            const std::set<const stk::mesh::Part*> &exclude)
00666 {
00667   Ioss::Region *region = mesh_data.m_output_region;
00668   region->begin_mode(Ioss::STATE_TRANSIENT);
00669 
00670   int out_step = region->add_state(time);
00671   internal_process_output_request(mesh_data, bulk, out_step,exclude);
00672 
00673   region->end_mode(Ioss::STATE_TRANSIENT);
00674 
00675   return out_step;
00676 }
00677 
00678 // ========================================================================
00679 void populate_bulk_data(stk::mesh::BulkData &bulk_data,
00680                         MeshData &mesh_data)
00681 {
00682   Ioss::Region *region = mesh_data.m_input_region;
00683 
00684   if (region) {
00685 
00686     bulk_data.modification_begin();
00687     process_mesh_bulk_data(region, bulk_data);
00688     bulk_data.modification_end();
00689 
00690   } else {
00691   std::cerr << "INTERNAL ERROR: Mesh Input Region pointer is NULL in populate_bulk_data.\n";
00692   std::exit(EXIT_FAILURE);
00693   }
00694 }
00695 
00696 // ========================================================================
00697 void process_mesh_bulk_data(Ioss::Region *region, stk::mesh::BulkData &bulk_data)
00698 {
00699     bool ints64bit = db_api_int_size(region) == 8;
00700     if (ints64bit) {
00701       int64_t zero = 0;
00702       process_elementblocks(*region, bulk_data, zero);
00703       process_nodeblocks(*region,    bulk_data);
00704       process_nodesets(*region,      bulk_data, zero);
00705       process_sidesets(*region,      bulk_data);
00706     } else {
00707       int zero = 0;
00708       process_elementblocks(*region, bulk_data, zero);
00709       process_nodeblocks(*region,    bulk_data);
00710       process_nodesets(*region,      bulk_data, zero);
00711       process_sidesets(*region,      bulk_data);
00712     }
00713 }
00714 
00715 namespace {
00716 // ========================================================================
00717 // Transfer transient field data from mesh file for io_entity to
00718 // the corresponding stk_mesh entities If there is a stk_mesh
00719 // field with the same name as the database field.
00720 // Assumes that mesh is positioned at the correct state for reading.
00721 void internal_process_input_request(Ioss::GroupingEntity *io_entity,
00722                                     stk::mesh::EntityRank entity_rank,
00723                                     stk::mesh::BulkData &bulk)
00724 {
00725   assert(io_entity != NULL);
00726   std::vector<stk::mesh::Entity*> entity_list;
00727   stk::io::get_entity_list(io_entity, entity_rank, bulk, entity_list);
00728 
00729   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00730 
00731   Ioss::NameList names;
00732   io_entity->field_describe(Ioss::Field::TRANSIENT, &names);
00733   for (Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00734     stk::mesh::FieldBase *field = fem_meta.get_field<stk::mesh::FieldBase>(*I);
00735     if (field) {
00736       stk::io::field_data_from_ioss(field, entity_list, io_entity, *I);
00737     }
00738   }
00739 }
00740 
00741 void input_nodeblock_fields(Ioss::Region &region, stk::mesh::BulkData &bulk)
00742 {
00743   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00744   const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
00745   assert(node_blocks.size() == 1);
00746 
00747   Ioss::NodeBlock *nb = node_blocks[0];
00748   internal_process_input_request(nb, fem_meta.node_rank(), bulk);
00749 }
00750 
00751 void input_elementblock_fields(Ioss::Region &region, stk::mesh::BulkData &bulk)
00752 {
00753   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00754   const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
00755   for(size_t i=0; i < elem_blocks.size(); i++) {
00756     if (stk::io::include_entity(elem_blocks[i])) {
00757       internal_process_input_request(elem_blocks[i], fem_meta.element_rank(), bulk);
00758     }
00759   }
00760 }
00761 
00762 void input_nodeset_fields(Ioss::Region &region, stk::mesh::BulkData &bulk)
00763 {
00764   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00765   const Ioss::NodeSetContainer& nodesets = region.get_nodesets();
00766   for(size_t i=0; i < nodesets.size(); i++) {
00767     if (stk::io::include_entity(nodesets[i])) {
00768       internal_process_input_request(nodesets[i], fem_meta.node_rank(), bulk);
00769     }
00770   }
00771 }
00772 
00773 void input_sideset_fields(Ioss::Region &region, stk::mesh::BulkData &bulk)
00774 {
00775   const stk::mesh::fem::FEMMetaData &fem_meta = stk::mesh::fem::FEMMetaData::get(bulk);
00776   if (fem_meta.spatial_dimension() <= fem_meta.side_rank())
00777     return;
00778 
00779   const Ioss::SideSetContainer& side_sets = region.get_sidesets();
00780   for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
00781       it != side_sets.end(); ++it) {
00782     Ioss::SideSet *entity = *it;
00783     if (stk::io::include_entity(entity)) {
00784       const Ioss::SideBlockContainer& blocks = entity->get_side_blocks();
00785       for(size_t i=0; i < blocks.size(); i++) {
00786         if (stk::io::include_entity(blocks[i])) {
00787           internal_process_input_request(blocks[i], fem_meta.side_rank(), bulk);
00788         }
00789       }
00790     }
00791   }
00792 }
00793 
00794 void define_input_nodeblock_fields(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00795 {
00796   const Ioss::NodeBlockContainer& node_blocks = region.get_node_blocks();
00797   assert(node_blocks.size() == 1);
00798 
00799   Ioss::NodeBlock *nb = node_blocks[0];
00800   stk::io::define_io_fields(nb, Ioss::Field::TRANSIENT,
00801                             fem_meta.universal_part(), fem_meta.node_rank());
00802 }
00803 
00804 void define_input_elementblock_fields(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00805 {
00806   const Ioss::ElementBlockContainer& elem_blocks = region.get_element_blocks();
00807   for(size_t i=0; i < elem_blocks.size(); i++) {
00808     if (stk::io::include_entity(elem_blocks[i])) {
00809       stk::mesh::Part* const part = fem_meta.get_part(elem_blocks[i]->name());
00810       assert(part != NULL);
00811       stk::io::define_io_fields(elem_blocks[i], Ioss::Field::TRANSIENT,
00812                                 *part, part_primary_entity_rank(*part));
00813     }
00814   }
00815 }
00816 
00817 void define_input_nodeset_fields(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00818 {
00819   const Ioss::NodeSetContainer& nodesets = region.get_nodesets();
00820   for(size_t i=0; i < nodesets.size(); i++) {
00821     if (stk::io::include_entity(nodesets[i])) {
00822       stk::mesh::Part* const part = fem_meta.get_part(nodesets[i]->name());
00823       assert(part != NULL);
00824       stk::io::define_io_fields(nodesets[i], Ioss::Field::TRANSIENT,
00825                                 *part, part_primary_entity_rank(*part));
00826     }
00827   }
00828 }
00829 
00830 void define_input_sideset_fields(Ioss::Region &region, stk::mesh::fem::FEMMetaData &fem_meta)
00831 {
00832   if (fem_meta.spatial_dimension() <= fem_meta.side_rank())
00833     return;
00834 
00835   const Ioss::SideSetContainer& side_sets = region.get_sidesets();
00836   for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
00837       it != side_sets.end(); ++it) {
00838     Ioss::SideSet *entity = *it;
00839     if (stk::io::include_entity(entity)) {
00840       const Ioss::SideBlockContainer& blocks = entity->get_side_blocks();
00841       for(size_t i=0; i < blocks.size(); i++) {
00842         if (stk::io::include_entity(blocks[i])) {
00843           stk::mesh::Part* const part = fem_meta.get_part(blocks[i]->name());
00844           assert(part != NULL);
00845           stk::io::define_io_fields(blocks[i], Ioss::Field::TRANSIENT,
00846                                     *part, part_primary_entity_rank(*part));
00847         }
00848       }
00849     }
00850   }
00851 }
00852 
00853 }
00854 
00855 // ========================================================================
00856 // Iterate over all Ioss entities in the input mesh database and
00857 // define a stk_field for all transient fields found.  The stk
00858 // field will have the same name as the field on the database.
00859 //
00860 // Note that all fields found on the database will have a
00861 // corresponding stk field defined.  If you want just a selected
00862 // subset of the defined fields, you will need to define the
00863 // fields manually.
00864 //
00865 // To populate the stk field with data from the database, call
00866 // process_input_request().
00867 void define_input_fields(MeshData &mesh_data, stk::mesh::fem::FEMMetaData &fem_meta)
00868 {
00869   Ioss::Region *region = mesh_data.m_input_region;
00870   if (region) {
00871   define_input_nodeblock_fields(*region, fem_meta);
00872   define_input_elementblock_fields(*region, fem_meta);
00873   define_input_nodeset_fields(*region, fem_meta);
00874   define_input_sideset_fields(*region, fem_meta);
00875   } else {
00876   std::cerr << "INTERNAL ERROR: Mesh Input Region pointer is NULL in process_input_request.\n";
00877   std::exit(EXIT_FAILURE);
00878   }
00879 }
00880 
00881 // ========================================================================
00882 // Iterate over all fields defined in the stk mesh data structure.
00883 // If the field has the io_attribute set, then define that field
00884 // on the corresponding io entity on the output mesh database.
00885 // The database field will have the same name as the stk field.
00886 //
00887 // To export the data to the database, call
00888 // process_output_request().
00889 
00890 void define_output_fields(const MeshData &mesh_data, const stk::mesh::fem::FEMMetaData &fem_meta,
00891                           bool add_all_fields)
00892 {
00893   Ioss::Region *region = mesh_data.m_output_region;
00894   if (region) {
00895   region->begin_mode(Ioss::STATE_DEFINE_TRANSIENT);
00896 
00897   // Special processing for nodeblock (all nodes in model)...
00898   stk::io::ioss_add_fields(fem_meta.universal_part(), fem_meta.node_rank(),
00899                              region->get_node_blocks()[0],
00900                              Ioss::Field::TRANSIENT, add_all_fields);
00901 
00902   const stk::mesh::PartVector & all_parts = fem_meta.get_parts();
00903   for ( stk::mesh::PartVector::const_iterator
00904             ip = all_parts.begin(); ip != all_parts.end(); ++ip ) {
00905 
00906     stk::mesh::Part * const part = *ip;
00907 
00908     // Check whether this part should be output to results database.
00909     if (stk::io::is_part_io_part(*part)) {
00910       // Get Ioss::GroupingEntity corresponding to this part...
00911       Ioss::GroupingEntity *entity = region->get_entity(part->name());
00912       if (entity != NULL) {
00913         stk::io::ioss_add_fields(*part, part_primary_entity_rank(*part),
00914                                    entity, Ioss::Field::TRANSIENT, add_all_fields);
00915       }
00916     }
00917   }
00918   region->end_mode(Ioss::STATE_DEFINE_TRANSIENT);
00919   } else {
00920   std::cerr << "INTERNAL ERROR: Mesh Input Region pointer is NULL in process_input_request.\n";
00921   std::exit(EXIT_FAILURE);
00922   }
00923 }
00924 // ========================================================================
00925 void process_input_request(MeshData &mesh_data, stk::mesh::BulkData &bulk, double time)
00926 {
00927   // Find the step on the database with time closest to the requested time...
00928   Ioss::Region *region = mesh_data.m_input_region;
00929   int step_count = region->get_property("state_count").get_int();
00930   double delta_min = 1.0e30;
00931   int    step_min  = 0;
00932   for (int istep = 0; istep < step_count; istep++) {
00933   double state_time = region->get_state_time(istep+1);
00934   double delta = state_time - time;
00935   if (delta < 0.0) delta = -delta;
00936   if (delta < delta_min) {
00937     delta_min = delta;
00938     step_min  = istep;
00939     if (delta == 0.0) break;
00940   }
00941   }
00942   // Exodus steps are 1-based;
00943   process_input_request(mesh_data, bulk, step_min+1);
00944 }
00945 
00946 void process_input_request(MeshData &mesh_data,
00947                            stk::mesh::BulkData &bulk,
00948                            int step)
00949 {
00950   if (step <= 0)
00951   return;
00952 
00953   Ioss::Region *region = mesh_data.m_input_region;
00954   if (region) {
00955   bulk.modification_begin();
00956 
00957   input_mesh_fields(region, bulk, step);
00958 
00959   bulk.modification_end();
00960 
00961   } else {
00962   std::cerr << "INTERNAL ERROR: Mesh Input Region pointer is NULL in process_input_request.\n";
00963   std::exit(EXIT_FAILURE);
00964   }
00965 }
00966 
00967 void input_mesh_fields(Ioss::Region *region, stk::mesh::BulkData &bulk,
00968                            double time)
00969 {
00970   // Find the step on the database with time closest to the requested time...
00971   int step_count = region->get_property("state_count").get_int();
00972   double delta_min = 1.0e30;
00973   int    step_min  = 0;
00974   for (int istep = 0; istep < step_count; istep++) {
00975   double state_time = region->get_state_time(istep+1);
00976   double delta = state_time - time;
00977   if (delta < 0.0) delta = -delta;
00978   if (delta < delta_min) {
00979     delta_min = delta;
00980     step_min  = istep;
00981     if (delta == 0.0) break;
00982   }
00983   }
00984   // Exodus steps are 1-based;
00985   input_mesh_fields(region, bulk, step_min+1);
00986 }
00987 
00988 void input_mesh_fields(Ioss::Region *region, stk::mesh::BulkData &bulk,
00989                            int step)
00990 {
00991   // Pick which time index to read into solution field.
00992   region->begin_state(step);
00993 
00994   input_nodeblock_fields(*region, bulk);
00995   input_elementblock_fields(*region, bulk);
00996   input_nodeset_fields(*region, bulk);
00997   input_sideset_fields(*region, bulk);
00998 
00999   region->end_state(step);
01000 }
01001 
01002 // ========================================================================
01003 template <typename INT>
01004 void get_element_block_sizes(MeshData &mesh_data,
01005                              std::vector<INT>& el_blocks)
01006 {
01007   Ioss::Region *io = mesh_data.m_input_region;
01008   const Ioss::ElementBlockContainer& elem_blocks = io->get_element_blocks();
01009   for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin(); it != elem_blocks.end(); ++it) {
01010     Ioss::ElementBlock *entity = *it;
01011     if (stk::io::include_entity(entity)) {
01012       el_blocks.push_back(entity->get_property("entity_count").get_int());
01013     }
01014   }
01015 }
01016 template void get_element_block_sizes(MeshData &mesh_data, std::vector<int>& el_blocks);
01017 template void get_element_block_sizes(MeshData &mesh_data, std::vector<int64_t>& el_blocks);
01018 } // namespace io
01019 } // namespace stk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines