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