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