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