Sierra Toolkit Version of the Day
IossBridge.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 <string.h>
00010 #include <iostream>
00011 #include <complex>
00012 
00013 #include <init/Ionit_Initializer.h>
00014 #include <Ioss_SubSystem.h>
00015 #include <Ioss_NullEntity.h>
00016 
00017 #include <stk_util/util/tokenize.hpp>
00018 #include <stk_io/IossBridge.hpp>
00019 
00020 #include <stk_util/parallel/Parallel.hpp>
00021 
00022 #include <stk_mesh/base/Types.hpp>
00023 #include <stk_mesh/base/Field.hpp>
00024 #include <stk_mesh/base/GetEntities.hpp>
00025 #include <stk_mesh/base/FieldData.hpp>
00026 #include <stk_mesh/base/MetaData.hpp>
00027 #include <stk_mesh/base/BulkData.hpp>
00028 #include <stk_mesh/fem/FEMHelpers.hpp>
00029 
00030 #include <Shards_BasicTopologies.hpp>
00031 
00032 namespace {
00033 
00034 stk::mesh::EntityRank get_entity_rank(const Ioss::GroupingEntity *entity,
00035                                       const stk::mesh::MetaData &meta)
00036 {
00037   switch (entity->type()) {
00038   case Ioss::NODEBLOCK:
00039     return stk::io::node_rank(meta);
00040 
00041   case Ioss::NODESET:
00042     return stk::io::node_rank(meta);
00043 
00044   case Ioss::ELEMENTBLOCK:
00045     return stk::io::element_rank(meta);
00046 
00047   case Ioss::SUPERELEMENT:
00048     return stk::io::element_rank(meta);
00049 
00050   case Ioss::SIDESET:
00051     {
00052       const Ioss::SideSet *sset = dynamic_cast<const Ioss::SideSet*>(entity);
00053       assert(sset != NULL);
00054       int my_rank = sset->max_parametric_dimension();
00055       if (my_rank == 2)
00056         return stk::io::face_rank(meta);
00057       if (my_rank == 1)
00058         return stk::io::edge_rank(meta);
00059       if (my_rank == 0)
00060         return stk::io::node_rank(meta);
00061       else
00062         return stk::mesh::InvalidEntityRank;
00063     }
00064 
00065   case Ioss::SIDEBLOCK:
00066     {
00067       const Ioss::SideBlock *sblk = dynamic_cast<const Ioss::SideBlock*>(entity);
00068       assert(sblk != NULL);
00069       int rank = sblk->topology()->parametric_dimension();
00070       if (rank == 2)
00071         return stk::io::face_rank(meta);
00072       if (rank == 1)
00073         return stk::io::edge_rank(meta);
00074       if (rank == 0)
00075         return stk::io::node_rank(meta);
00076       else
00077         return stk::mesh::InvalidEntityRank;
00078     }
00079   default:
00080     return stk::mesh::InvalidEntityRank;
00081   }
00082 }
00083 
00084 const CellTopologyData *map_ioss_to_topology( const std::string &element_type ,
00085                                               const int node_count)
00086 {
00089 
00090   const CellTopologyData* celltopo = NULL ;
00091 
00092   const char *etype = element_type.c_str();
00093   if ( 0 == strncasecmp( "circle" , etype , 6 ) ) {
00094     celltopo = shards::getCellTopologyData< shards::Particle >();
00095   }
00096   else if ( 0 == strncasecmp( "sphere" , etype , 6) ) {
00097     celltopo = shards::getCellTopologyData< shards::Particle >();
00098   }
00099   // bar, beam, truss, rod...
00100   else if ( 0 == strncasecmp( "bar" , etype , 3 ) ) {
00101     if ( node_count == 2 ) {
00102       celltopo = shards::getCellTopologyData< shards::Beam<2> >();
00103     }
00104     else if ( node_count == 3 ) {
00105       celltopo = shards::getCellTopologyData< shards::Beam<3> >();
00106     }
00107   }
00108   else if ( 0 == strncasecmp( "shellline2d" , etype , 11 ) ) {
00109     if ( node_count == 2) {
00110       celltopo = shards::getCellTopologyData< shards::ShellLine<2> >();
00111     }
00112     else if ( node_count == 3) {
00113       celltopo = shards::getCellTopologyData< shards::ShellLine<3> >();
00114     }
00115   } else if ( 0 == strncasecmp( "shell" , etype , 5 ) ) {
00116     // shell4, shell8, shell9
00117     if ( node_count == 4 ) {
00118       celltopo = shards::getCellTopologyData< shards::ShellQuadrilateral<4> >();
00119     }
00120     else if ( node_count == 8 ) {
00121       celltopo = shards::getCellTopologyData< shards::ShellQuadrilateral<8> >();
00122     }
00123     else if ( node_count == 9 ) {
00124       celltopo = shards::getCellTopologyData< shards::ShellQuadrilateral<9> >();
00125     }
00126   }
00127   else if ( 0 == strncasecmp( "quad" , etype , 3 ) ) {
00128     // The 2D types would be quad4, quad8, and quad9.
00129     // The 3D types would be quad faces of a hex... quadface4,
00130     // quadface8, quadface9.
00131     if ( node_count == 4 ) {
00132       celltopo = shards::getCellTopologyData< shards::Quadrilateral<4> >();
00133     }
00134     else if ( node_count == 8 ) {
00135       celltopo = shards::getCellTopologyData< shards::Quadrilateral<8> >();
00136     }
00137     else if ( node_count == 9 ) {
00138       celltopo = shards::getCellTopologyData< shards::Quadrilateral<9> >();
00139     }
00140   }
00141   else if ( 0 == strncasecmp( "trishell" , etype , 8 ) ) {
00142     if ( node_count == 3 ) {
00143       celltopo = shards::getCellTopologyData< shards::ShellTriangle<3> >();
00144     }
00145     else if ( node_count == 6 ) {
00146       celltopo = shards::getCellTopologyData< shards::ShellTriangle<6> >();
00147     }
00148   }
00149 
00150   else if (0 == strncasecmp("triface", etype, 7) ||
00151            0 == strncasecmp("tri",     etype, 3)) {
00152     if ( node_count == 3 ) {
00153       celltopo = shards::getCellTopologyData< shards::Triangle<3> >();
00154     }
00155     else if ( node_count == 4 ) {
00156       celltopo = shards::getCellTopologyData< shards::Triangle<4> >();
00157     }
00158     else if ( node_count == 6 ) {
00159       celltopo = shards::getCellTopologyData< shards::Triangle<6> >();
00160     }
00161   }
00162 
00163   else if ( 0 == strncasecmp( "pyramid" , etype , 7 ) ) {
00164     if ( node_count == 5 ) {
00165       celltopo = shards::getCellTopologyData< shards::Pyramid<5> >();
00166     }
00167     else if ( node_count == 13 ) {
00168       celltopo = shards::getCellTopologyData< shards::Pyramid<13> >();
00169     }
00170     else if ( node_count == 14 ) {
00171       celltopo = shards::getCellTopologyData< shards::Pyramid<14> >();
00172     }
00173   }
00174 
00176   else if ( 0 == strncasecmp( "tetra" , etype , 5 ) ) {
00177     if ( node_count == 4 ) {
00178       celltopo = shards::getCellTopologyData< shards::Tetrahedron<4> >();
00179     }
00180     else if ( node_count ==  8 ) {
00181       celltopo = shards::getCellTopologyData< shards::Tetrahedron<8> >();
00182     }
00183     else if ( node_count == 10 ) {
00184       celltopo = shards::getCellTopologyData< shards::Tetrahedron<10> >();
00185     }
00186   }
00187   else if ( 0 == strncasecmp( "wedge" , etype , 5 ) ) {
00188     if ( node_count == 6 ) {
00189       celltopo = shards::getCellTopologyData< shards::Wedge<6> >();
00190     }
00191     else if ( node_count == 15 ) {
00192       celltopo = shards::getCellTopologyData< shards::Wedge<15> >();
00193     }
00194     else if ( node_count == 18 ) {
00195       celltopo = shards::getCellTopologyData< shards::Wedge<18> >();
00196     }
00197   }
00198   else if ( 0 == strncasecmp( "hex" , etype , 3 ) ) {
00199     if ( node_count == 8 ) {
00200       celltopo = shards::getCellTopologyData< shards::Hexahedron<8> >();
00201     }
00202     else if ( node_count == 20 ) {
00203       celltopo = shards::getCellTopologyData< shards::Hexahedron<20> >();
00204     }
00205     else if ( node_count == 27 ) {
00206       celltopo = shards::getCellTopologyData< shards::Hexahedron<27> >();
00207     }
00208   }
00209 
00210   else if (0 == strncasecmp("edge", etype, 4)) {
00211     if ( node_count == 2) {
00212       // edge2, edge2d2, edge3d2
00213       celltopo = shards::getCellTopologyData< shards::Line<2> >();
00214     }
00215     else if ( node_count == 3) {
00216       // edge3, edge2d3, edge3d3
00217       celltopo = shards::getCellTopologyData< shards::Line<3> >();
00218     }
00219   }
00220 
00221   else if (0 == strncasecmp("node", etype, 4)) {
00222     celltopo = shards::getCellTopologyData< shards::Node >();
00223   }
00224 
00225   if ( NULL == celltopo ) {
00226     std::ostringstream oss;
00227     oss << "ERROR, unsupported topology name = '" << element_type
00228         << "' , node_count = " << node_count;
00229     throw std::runtime_error(oss.str());
00230   }
00231 
00232   return celltopo;
00233 }
00234 
00235 template <typename T>
00236 const stk::mesh::FieldBase *declare_ioss_field_internal(stk::mesh::MetaData &meta,
00237                                                         stk::mesh::EntityRank type,
00238                                                         stk::mesh::Part &part,
00239                                                         const Ioss::Field &io_field,
00240                                                         bool use_cartesian_for_scalar, T /*dummy*/)
00241 {
00242   stk::mesh::FieldBase *field_ptr = NULL;
00243   std::string field_type = io_field.transformed_storage()->name();
00244   std::string name = io_field.get_name();
00245   size_t num_components = io_field.transformed_storage()->component_count();
00246 
00247   if (field_type == "scalar" || num_components == 1) {
00248     if (!use_cartesian_for_scalar) {
00249       stk::mesh::Field<double> & field = meta.declare_field<stk::mesh::Field<double> >(name);
00250       stk::mesh::put_field(field, type, part);
00251       field_ptr = &field;
00252     } else {
00253       stk::mesh::Field<double, stk::mesh::Cartesian> & field =
00254         meta.declare_field<stk::mesh::Field<double, stk::mesh::Cartesian> >(name);
00255       stk::mesh::put_field(field, type, part, 1);
00256       field_ptr = &field;
00257     }
00258   }
00259   else if (field_type == "vector_2d") {
00260     stk::mesh::Field<double, stk::mesh::Cartesian> & field =
00261       meta.declare_field<stk::mesh::Field<double, stk::mesh::Cartesian> >(name);
00262     stk::mesh::put_field(field, type, part, 2);
00263     field_ptr = &field;
00264   }
00265   else if (field_type == "vector_3d") {
00266     stk::mesh::Field<double, stk::mesh::Cartesian> & field =
00267       meta.declare_field<stk::mesh::Field<double,
00268       stk::mesh::Cartesian> >(name);
00269     stk::mesh::put_field(field, type, part, 3);
00270     field_ptr = &field;
00271   }
00272   else if (field_type == "sym_tensor_33") {
00273     stk::mesh::Field<double, stk::mesh::SymmetricTensor> & field =
00274       meta.declare_field<stk::mesh::Field<double,
00275       stk::mesh::SymmetricTensor> >(name);
00276     stk::mesh::put_field(field, type, part, 6);
00277     field_ptr = &field;
00278   }
00279   else if (field_type == "full_tensor_36") {
00280     stk::mesh::Field<double, stk::mesh::FullTensor> & field =
00281       meta.declare_field<stk::mesh::Field<double,
00282       stk::mesh::FullTensor> >(name);
00283     stk::mesh::put_field(field, type, part, 9);
00284     field_ptr = &field;
00285   }
00286   else {
00287     // Just create a field with the correct number of components...
00288     stk::mesh::Field<double,shards::ArrayDimension> & field =
00289       meta.declare_field<stk::mesh::Field<double,shards::ArrayDimension> >(name);
00290     stk::mesh::put_field(field, type, part, num_components);
00291     field_ptr = &field;
00292   }
00293 
00294   if (field_ptr != NULL) {
00295     stk::io::set_field_role(*field_ptr, io_field.get_role());
00296   }
00297   return field_ptr;
00298 }
00299 
00300 const stk::mesh::FieldBase *declare_ioss_field(stk::mesh::MetaData &meta,
00301                                                stk::mesh::EntityRank type,
00302                                                stk::mesh::Part &part,
00303                                                const Ioss::Field &io_field,
00304                                                bool use_cartesian_for_scalar)
00305 {
00306   const stk::mesh::FieldBase *field_ptr = NULL;
00307   if (io_field.get_type() == Ioss::Field::INTEGER) {
00308     field_ptr = declare_ioss_field_internal(meta, type, part, io_field, use_cartesian_for_scalar, (int)1);
00309   } else if (io_field.get_type() == Ioss::Field::REAL) {
00310     field_ptr = declare_ioss_field_internal(meta, type, part, io_field, use_cartesian_for_scalar, (double)1.0);
00311   } else if (io_field.get_type() == Ioss::Field::COMPLEX) {
00312     field_ptr = declare_ioss_field_internal(meta, type, part, io_field, use_cartesian_for_scalar, std::complex<double>(0.0,0.0));
00313   } else {
00314     std::ostringstream errmsg;
00315     errmsg << "ERROR: Unrecognized field type for IO field '"
00316            << io_field.get_name() << "'.";
00317     throw std::runtime_error(errmsg.str());
00318   }
00319   return field_ptr;
00320 }
00321 
00322 template <typename T>
00323 void internal_field_data_from_ioss(const Ioss::Field &io_field,
00324                                    const stk::mesh::FieldBase *field,
00325                                    std::vector<stk::mesh::Entity*> &entities,
00326                                    Ioss::GroupingEntity *io_entity,
00327                                    T /*dummy */)
00328 {
00329   size_t field_component_count = io_field.transformed_storage()->component_count();
00330 
00331   std::vector<T> io_field_data;
00332   size_t io_entity_count = io_entity->get_field_data(io_field.get_name(), io_field_data);
00333   assert(io_field_data.size() == entities.size() * field_component_count);
00334 
00335   size_t entity_count = entities.size();
00336 
00337   if (io_entity_count != entity_count) {
00338     std::ostringstream errmsg;
00339     errmsg << "ERROR: Field count mismatch for IO field '"
00340            << io_field.get_name() << "'. The IO system has " << io_entity_count
00341            << " entries, but the stk:mesh system has " << entity_count
00342            << " entries. The two counts must match.";
00343     throw std::runtime_error(errmsg.str());
00344   }
00345 
00346   for (size_t i=0; i < entity_count; ++i) {
00349     if (entities[i] != NULL) {
00350       T *fld_data = (T*)stk::mesh::field_data(*field, *entities[i]);
00351       assert(fld_data != NULL);
00352       for(size_t j=0; j<field_component_count; ++j) {
00353         fld_data[j] = io_field_data[i*field_component_count+j];
00354       }
00355     }
00356   }
00357 }
00358 
00359 template <typename T>
00360 void internal_field_data_to_ioss(const Ioss::Field &io_field,
00361                                  const stk::mesh::FieldBase *field,
00362                                  std::vector<stk::mesh::Entity*> &entities,
00363                                  Ioss::GroupingEntity *io_entity,
00364                                  T /*dummy */)
00365 {
00366   size_t field_component_count = io_field.transformed_storage()->component_count();
00367   size_t entity_count = entities.size();
00368 
00369   std::vector<T> io_field_data(entity_count*field_component_count);
00370 
00371   for (size_t i=0; i < entity_count; ++i) {
00374     if (entities[i] != NULL) {
00375       T *fld_data = (T*)stk::mesh::field_data(*field, *entities[i]);
00376       assert(fld_data != NULL);
00377       for(size_t j=0; j<field_component_count; ++j) {
00378         io_field_data[i*field_component_count+j] = fld_data[j];
00379       }
00380     } else {
00381       for(size_t j=0; j<field_component_count; ++j) {
00382         io_field_data[i*field_component_count+j] = 0;
00383       }
00384     }
00385   }
00386 
00387   size_t io_entity_count = io_entity->put_field_data(io_field.get_name(), io_field_data);
00388   assert(io_field_data.size() == entities.size() * field_component_count);
00389 
00390   if (io_entity_count != entity_count) {
00391     std::ostringstream errmsg;
00392     errmsg << "ERROR: Field count mismatch for IO field '"
00393            << io_field.get_name() << "'. The IO system has " << io_entity_count
00394            << " entries, but the stk:mesh system has " << entity_count
00395            << " entries. The two counts must match.";
00396     throw std::runtime_error(errmsg.str());
00397   }
00398 }
00399 
00400 }//namespace <empty>
00401 
00402 namespace stk {
00403 namespace io {
00404 
00405 size_t db_api_int_size(const Ioss::GroupingEntity *entity)
00406 {
00407   return entity->get_database()->int_byte_size_api();
00408 }
00409 
00410 bool invalid_rank(stk::mesh::EntityRank rank)
00411 {
00412   return rank == mesh::InvalidEntityRank;
00413 }
00414 
00415 stk::mesh::EntityRank part_primary_entity_rank(const stk::mesh::Part &part)
00416 {
00417   if (mesh::MetaData::get(part).universal_part() == part) {
00418     return stk::mesh::fem::FEMMetaData::NODE_RANK;
00419   }
00420   else {
00421     return part.primary_entity_rank();
00422   }
00423 }
00424 
00425 stk::mesh::EntityRank element_rank(const stk::mesh::MetaData &meta)
00426 {
00427   return stk::mesh::fem::FEMMetaData::get(meta).element_rank();
00428 }
00429 
00430 stk::mesh::EntityRank side_rank(const stk::mesh::MetaData &meta)
00431 {
00432   return stk::mesh::fem::FEMMetaData::get(meta).side_rank();
00433 }
00434 
00435 stk::mesh::EntityRank face_rank(const stk::mesh::MetaData &meta)
00436 {
00437   return stk::mesh::fem::FEMMetaData::get(meta).face_rank();
00438 }
00439 
00440 stk::mesh::EntityRank edge_rank(const stk::mesh::MetaData &meta)
00441 {
00442   return stk::mesh::fem::FEMMetaData::get(meta).edge_rank();
00443 }
00444 
00445 stk::mesh::EntityRank node_rank(const stk::mesh::MetaData& meta)
00446 {
00447   return stk::mesh::fem::FEMMetaData::NODE_RANK;
00448 }
00449 
00450 void set_cell_topology(stk::mesh::Part &part, const CellTopologyData * const cell_topology)
00451 {
00452   stk::mesh::fem::set_cell_topology(part, cell_topology);
00453 }
00454 
00455 const CellTopologyData *get_cell_topology(const stk::mesh::Part &part)
00456 {
00457   return stk::mesh::fem::FEMMetaData::get(part).get_cell_topology(part).getCellTopologyData();
00458 }
00459 
00460 void initialize_spatial_dimension(stk::mesh::fem::FEMMetaData & fem_meta, size_t spatial_dimension,
00461                                   const std::vector<std::string> &entity_rank_names)
00462 {
00463   if (!fem_meta.is_FEM_initialized() ) {
00464     fem_meta.FEM_initialize(spatial_dimension, entity_rank_names);
00465   }
00466 }
00467 
00468 void initialize_spatial_dimension(stk::mesh::MetaData &meta, size_t spatial_dimension,
00469                                   const std::vector<std::string> &entity_rank_names)
00470 {
00471   stk::mesh::fem::FEMMetaData & fem_meta = stk::mesh::fem::FEMMetaData::get(meta);
00472   initialize_spatial_dimension(fem_meta, spatial_dimension, entity_rank_names);
00473 }
00474 
00475 void get_io_field_type(const stk::mesh::FieldBase *field,
00476                        const stk::mesh::FieldRestriction &res,
00477                        std::pair<std::string, Ioss::Field::BasicType> *result)
00478 {
00479   static const std::string invalid("invalid");
00480   static const std::string scalar("scalar");
00481   static const std::string vector_2d("vector_2d");
00482   static const std::string vector_3d("vector_3d");
00483   static const std::string quaternion_2d("quaternion_2d");
00484   static const std::string quaternion_3d("quaternion_3d");
00485   static const std::string full_tensor_36("full_tensor_36");
00486   static const std::string full_tensor_32("full_tensor_32");
00487   static const std::string full_tensor_22("full_tensor_22");
00488   static const std::string full_tensor_16("full_tensor_16");
00489   static const std::string full_tensor_12("full_tensor_12");
00490   static const std::string sym_tensor_33("sym_tensor_33");
00491   static const std::string sym_tensor_31("sym_tensor_31");
00492   static const std::string sym_tensor_21("sym_tensor_21");
00493   static const std::string sym_tensor_13("sym_tensor_13");
00494   static const std::string sym_tensor_11("sym_tensor_11");
00495   static const std::string sym_tensor_10("sym_tensor_10");
00496   static const std::string asym_tensor_03("asym_tensor_03");
00497   static const std::string asym_tensor_02("asym_tensor_02");
00498   static const std::string asym_tensor_01("asym_tensor_01");
00499   static const std::string matrix_22("matrix_22");
00500   static const std::string matrix_33("matrix_33");
00501 
00502   const unsigned rank = field->rank();
00503   const shards::ArrayDimTag * const * const tags = field->dimension_tags();
00504 
00505   result->second = Ioss::Field::INVALID;
00506 
00507   if ( field->type_is<double>() ) {
00508   result->second = Ioss::Field::REAL;
00509   }
00510   else if ( field->type_is<int>() ) {
00511   result->second = Ioss::Field::INTEGER;
00512   }
00513 
00514   if ( 0 == rank ) {
00515   result->first = scalar ;
00516   }
00517   else if ( 1 == rank ) {
00518   size_t num_comp = res.stride(0);
00519   if ( tags[0] == & stk::mesh::Cartesian::tag() && 1 == num_comp ) {
00520     result->first = scalar ;
00521   }
00522   else if ( tags[0] == & stk::mesh::Cartesian::tag() && 2 == num_comp ) {
00523     result->first = vector_2d ;
00524   }
00525   else if ( tags[0] == & stk::mesh::Cartesian::tag() && 3 == num_comp ) {
00526     result->first = vector_3d ;
00527   }
00528   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 9 == num_comp ) {
00529     result->first = full_tensor_36 ;
00530   }
00531   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 5 == num_comp ) {
00532     result->first = full_tensor_32 ;
00533   }
00534   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 4 == num_comp ) {
00535     result->first = full_tensor_22 ;
00536   }
00537   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 3 == num_comp ) {
00538     result->first = full_tensor_12 ;
00539   }
00540   else if ( tags[0] == & stk::mesh::SymmetricTensor::tag() && 6 == num_comp ) {
00541     result->first = sym_tensor_33 ;
00542   }
00543   else if ( tags[0] == & stk::mesh::SymmetricTensor::tag() && 4 == num_comp ) {
00544     result->first = sym_tensor_31 ;
00545   }
00546   else if ( tags[0] == & stk::mesh::SymmetricTensor::tag() && 3 == num_comp ) {
00547     result->first = sym_tensor_21 ;
00548   }
00549   }
00550 
00551   if ( result->first.empty() ) {
00552   size_t num_comp = res.stride(rank-1);
00553   std::ostringstream tmp ;
00554   tmp << "Real[" << num_comp << "]" ;
00555   result->first = tmp.str();
00556   }
00557 }
00558 
00559 //----------------------------------------------------------------------
00560 const Ioss::GroupingEntity *get_associated_ioss_entity(const mesh::Part &part)
00561 {
00562   const Ioss::GroupingEntity *entity = part.attribute<Ioss::GroupingEntity>();
00563   if (!entity || entity->type() == Ioss::INVALID_TYPE) {
00564   return NULL;
00565   } else {
00566   return entity;
00567   }
00568 }
00569 
00570 void put_io_part_attribute(mesh::Part & part, Ioss::GroupingEntity *entity)
00571 {
00572   if (part.attribute<Ioss::GroupingEntity>() != NULL) {
00573   std::string msg = "stk::io::put_io_part_attribute( ";
00574   msg += part.name();
00575   msg += " ) FAILED:";
00576   msg += " io_part_attribute is already defined";
00577   throw std::runtime_error( msg );
00578   }
00579 
00580   mesh::MetaData & meta = mesh::MetaData::get(part);
00581   if (entity) {
00582   meta.declare_attribute_no_delete(part, entity);
00583   } else {
00584   Ioss::GroupingEntity *attr = new Ioss::NullEntity();
00585   meta.declare_attribute_with_delete(part, attr);
00586   }
00587 }
00588 
00589 void remove_io_part_attribute(mesh::Part & part)
00590 {
00591   const Ioss::GroupingEntity *entity = part.attribute<Ioss::GroupingEntity>();
00592   if (entity == NULL) {
00593   std::string msg = "stk::io::remove_io_part_attribute( ";
00594   msg += part.name();
00595   msg += " ) FAILED:";
00596   msg += " io_part_attribute is not defined on this part";
00597   throw std::runtime_error( msg );
00598   } else {
00599   mesh::MetaData & meta = mesh::MetaData::get(part);
00600   bool success = meta.remove_attribute(part, entity);
00601   if (!success) {
00602     std::string msg = "stk::io::remove_io_part_attribute( ";
00603     msg += part.name();
00604     msg += " ) FAILED:";
00605     msg += " meta.remove_attribute(..) returned failure.";
00606     throw std::runtime_error( msg );
00607   }
00608 
00609   if (entity->type() == Ioss::INVALID_TYPE) {
00610     delete entity;
00611   }
00612 
00613   }
00614 }
00615 
00620 bool is_valid_part_field(const stk::mesh::FieldBase *field,
00621                          const stk::mesh::EntityRank part_type,
00622                          const stk::mesh::Part &part,
00623                          const stk::mesh::Part &universal,
00624                          const Ioss::Field::RoleType filter_role,
00625                          bool add_all)
00626 {
00627   const Ioss::Field::RoleType *role = stk::io::get_field_role(*field);
00628 
00629   if (!add_all && role == NULL) {
00630   return false;
00631   }
00632 
00633   if (role != NULL && *role != filter_role)
00634   return false;
00635 
00636   const stk::mesh::FieldBase::Restriction &res = field->restriction(part_type, part);
00637   if (res.dimension() > 0) {
00638   // The field exists on the current 'part'.  Now check (for
00639   // node types only) whether the 'part' is *either* the
00640   // universal_part() *or* the field *doesn't* exist on the
00641   // universal part...
00642   // Note that for "node" type parts, the IO database has a part
00643   // (nodeblock) that corresponds to the universal_part(), so
00644   // fields defined on all nodes are output on the nodeblock and
00645   // fields defined on only a subset of the parts should be
00646   // output on that subset which maps to a nodeset.  For other
00647   // part types, there is no IO entity that corresponds to a
00648   // universal part, so fields must be output on the part they
00649   // exist on.  There may be a problem if we start using element
00650   // sets ..., but wait until we get to that point; current code
00651   // works with current entity set.
00652   if (part_type != node_rank(mesh::MetaData::get(part)) || part == universal) {
00653     return true;
00654   }
00655 
00656   const stk::mesh::FieldBase::Restriction &res_universe = field->restriction(part_type, universal);
00657   if (res_universe.dimension() <= 0) {
00658     // Field exists on current part, but not on the universal
00659     // set (and this part is not the universal part)
00660     return true;
00661   }
00662   }
00663   return false;
00664 }
00665 
00669 // ========================================================================
00670 const CellTopologyData *map_topology_ioss_to_cell(const Ioss::ElementTopology *topology)
00671 {
00678 
00679   std::string name         = topology->name();
00680   int io_nodes_per_element = topology->number_nodes();
00681 
00682   const CellTopologyData *cell_topology = map_ioss_to_topology(name, io_nodes_per_element);
00683 
00684   return cell_topology;
00685 }
00686 
00687 std::string map_topology_cell_to_ioss( const CellTopologyData *cell_top,
00688                                        int spatial_dimension)
00689 {
00690   std::string extype = "unknown";
00691 
00692   if (cell_top == NULL)
00693   return extype;
00694 
00695   if(strcmp(cell_top->name, "super") == 0) {
00696     std::stringstream oss;
00697     oss << "super" << cell_top->node_count;
00698     return oss.str();
00699   }
00700   else if(strncasecmp(cell_top->name, "super", 5) == 0) {
00701     return cell_top->name;
00702   }
00703 
00704   switch( cell_top->key ) {
00705   case shards::Node::key : extype.assign( "node" ); break ;
00706   case shards::Particle::key :
00707   if (spatial_dimension == 2) extype = "circle1";
00708   else                        extype = "sphere1";
00709   break ;
00710 
00711   case shards::Beam<2>::key :
00712   extype = "bar2";
00713   break ;
00714 
00715   case shards::Beam<3>::key :
00716   extype = "bar3";
00717   break ;
00718 
00719   case shards::Line<2>::key :
00720   if (spatial_dimension == 2) extype = "edge2d2";
00721   else                        extype = "edge2";
00722   break ;
00723   case shards::Line<3>::key :
00724   if (spatial_dimension == 2) extype = "edge2d3";
00725   else                        extype = "edge3";
00726   break ;
00727 
00728   case shards::ShellLine<2>::key : extype.assign( "shellline2d2" ); break ;
00729   case shards::ShellLine<3>::key : extype.assign( "shellline2d3" ); break ;
00730 
00731   case shards::Triangle<3>::key :
00732   extype = "tri3";
00733   break ;
00734   case shards::Triangle<4>::key :
00735     extype = "tri4";
00736   break ;
00737   case shards::Triangle<6>::key :
00738     extype = "tri6";
00739   break ;
00740 
00741   case shards::ShellTriangle<3>::key : extype.assign( "trishell3" ); break ;
00742   case shards::ShellTriangle<6>::key : extype.assign( "trishell6" ); break ;
00743 
00744   case shards::Quadrilateral<4>::key :
00745     extype = "quad4";
00746   break ;
00747   case shards::Quadrilateral<8>::key :
00748     extype = "quad8";
00749   break ;
00750   case shards::Quadrilateral<9>::key :
00751     extype = "quad9";
00752   break ;
00753 
00754   case shards::ShellQuadrilateral<4>::key : extype.assign( "shell4" ); break ;
00755   case shards::ShellQuadrilateral<8>::key : extype.assign( "shell8" ); break ;
00756   case shards::ShellQuadrilateral<9>::key : extype.assign( "shell9" ); break ;
00757 
00758   case shards::Tetrahedron< 4>::key : extype.assign( "tetra4" ); break ;
00759   case shards::Tetrahedron<8>::key : extype.assign( "tetra8" ); break ;
00760   case shards::Tetrahedron<10>::key : extype.assign( "tetra10" ); break ;
00761 
00762   case shards::Pyramid< 5>::key : extype.assign( "pyramid5" ); break ;
00763   case shards::Pyramid<13>::key : extype.assign( "pyramid13" ); break ;
00764   case shards::Pyramid<14>::key : extype.assign( "pyramid14" ); break ;
00765 
00766   case shards::Wedge< 6>::key : extype.assign( "wedge6" ); break ;
00767   case shards::Wedge<15>::key : extype.assign( "wedge15" ); break ;
00768   case shards::Wedge<18>::key : extype.assign( "wedge18" ); break ;
00769 
00770   case shards::Hexahedron< 8>::key : extype.assign( "hex8" ); break ;
00771   case shards::Hexahedron<20>::key : extype.assign( "hex20" ); break ;
00772   case shards::Hexahedron<27>::key : extype.assign( "hex27" ); break ;
00773 
00774   default:
00775   std::ostringstream oss;
00776   oss << "stk::io::map_topology_to_ioss( '" << cell_top->name
00777       << "' ) ERROR unmapped topology" << std::endl ;
00778   throw std::runtime_error(oss.str());
00779   }
00780 
00781   return extype ;
00782 }
00783 
00784 void internal_part_processing(Ioss::GroupingEntity *entity, stk::mesh::fem::FEMMetaData &meta)
00785 {
00786   internal_part_processing(entity, meta.get_meta_data(meta));
00787 }
00788 
00789 void internal_part_processing(Ioss::EntityBlock *entity, stk::mesh::fem::FEMMetaData &meta)
00790 {
00791   internal_part_processing(entity, meta.get_meta_data(meta));
00792 }
00793 
00794 void internal_part_processing(Ioss::GroupingEntity *entity, stk::mesh::MetaData &meta)
00795 {
00796   if (include_entity(entity)) {
00797     stk::mesh::fem::FEMMetaData * fem_meta = const_cast<stk::mesh::fem::FEMMetaData *>(meta.get_attribute<stk::mesh::fem::FEMMetaData>());
00798   mesh::EntityRank type = get_entity_rank(entity, meta);
00799     if (fem_meta) {
00800     stk::mesh::Part & part = fem_meta->declare_part(entity->name(), type);
00801     stk::io::put_io_part_attribute(part, entity);
00802     } else {
00803     stk::mesh::Part & part = meta.declare_part(entity->name(), type);
00804     stk::io::put_io_part_attribute(part, entity);
00805     }
00806   }
00807 }
00808 
00809 void internal_part_processing(Ioss::EntityBlock *entity, stk::mesh::MetaData &meta)
00810 {
00811   if (include_entity(entity)) {
00812   mesh::EntityRank type = get_entity_rank(entity, meta);
00813     stk::mesh::fem::FEMMetaData * fem_meta = const_cast<stk::mesh::fem::FEMMetaData *>(meta.get_attribute<stk::mesh::fem::FEMMetaData>());
00814     //const stk::mesh::fem::FEMMetaData * fem_meta = meta.get_attribute<stk::mesh::fem::FEMMetaData>();
00815     stk::mesh::Part * part = NULL;
00816     if( fem_meta )
00817       part = &fem_meta->declare_part(entity->name(), type);
00818     else
00819       part = &meta.declare_part(entity->name(), type);
00820     stk::io::put_io_part_attribute(*part, entity);
00821 
00822     const Ioss::ElementTopology *topology = entity->topology();
00823     // Check spatial dimension of the element topology here so we
00824     // can issue a more meaningful error message.  If the
00825     // dimension is bad and we continue to the following calls,
00826     // there is an exception and we get unintelligible (to the
00827     // user) error messages.  Could also do a catch...
00828 
00829     if (entity->type() == Ioss::ELEMENTBLOCK) {
00830       assert(topology != NULL);
00831       if (fem_meta && (topology->spatial_dimension() < (int)fem_meta->spatial_dimension())) {
00832   // NOTE: The comparison is '<' and not '!=' since a 2D mesh
00833   // can contain a "3d" element -- a Beam is both a 2D and
00834   // 3D element...
00835 
00836   std::ostringstream msg ;
00837   msg << "\n\nERROR: Element Block " << entity->name()
00838             << " contains " << topology->name() << " elements with spatial dimension "
00839             << topology->spatial_dimension()
00840             << "\n       which does not match the spatial dimension of the model which is "
00841             << fem_meta->spatial_dimension() << "\n\n";
00842       throw std::runtime_error( msg.str() );
00843       }
00844     }
00845 
00846     const CellTopologyData * const cell_topology = map_topology_ioss_to_cell(topology);
00847     // \todo IMPLEMENT Determine whether application can work
00848     // with this topology type... Perhaps map_topology_ioss_to_cell only
00849     // returns a valid topology if the application has registered
00850     // that it can handle that specific topology.
00851     
00852     if (cell_topology != NULL) {
00853       if( fem_meta ) {
00854   const stk::mesh::fem::CellTopology cell_topo(cell_topology);
00855   stk::mesh::fem::set_cell_topology(*part, cell_topo);
00856       }
00857       stk::io::set_cell_topology(*part, cell_topology);
00858     } else {
00859       // \todo IMPLEMENT handle cell_topolgy mapping error...
00860     }
00861     stk::io::define_io_fields(entity, Ioss::Field::ATTRIBUTE, *part, type);
00862   }
00863 }
00864 
00865 //----------------------------------------------------------------------
00870 void ioss_add_fields(const stk::mesh::Part &part,
00871                      const stk::mesh::EntityRank part_type,
00872                      Ioss::GroupingEntity *entity,
00873                      const Ioss::Field::RoleType filter_role,
00874                      const bool add_all)
00875 {
00876   const stk::mesh::MetaData & meta = mesh::MetaData::get(part);
00877   const stk::mesh::Part &universal = meta.universal_part();
00878 
00879   const std::vector<mesh::FieldBase*> &fields = meta.get_fields();
00880 
00881   std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
00882   while (I != fields.end()) {
00883   const stk::mesh::FieldBase *f = *I; ++I;
00884   if (stk::io::is_valid_part_field(f, part_type, part, universal, filter_role, add_all)) {
00885     const stk::mesh::FieldBase::Restriction &res = f->restriction(part_type, part);
00886     std::pair<std::string, Ioss::Field::BasicType> field_type;
00887     get_io_field_type(f, res, &field_type);
00888     if (field_type.second != Ioss::Field::INVALID) {
00889       size_t entity_size = entity->get_property("entity_count").get_int();
00890       const std::string& name = f->name();
00891       entity->field_add(Ioss::Field(name, field_type.second, field_type.first,
00892                                       filter_role, entity_size));
00893     }
00894   }
00895   }
00896 }
00897 
00907 void define_io_fields(Ioss::GroupingEntity *entity,
00908                       Ioss::Field::RoleType role,
00909                       stk::mesh::Part &part,
00910                       stk::mesh::EntityRank part_type)
00911 {
00912   stk::mesh::MetaData &meta = mesh::MetaData::get(part);
00913 
00914   bool use_cartesian_for_scalar = false;
00915   if (role == Ioss::Field::ATTRIBUTE)
00916   use_cartesian_for_scalar = true;
00917 
00918   Ioss::NameList names;
00919   entity->field_describe(role, &names);
00920 
00921   for (Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00922   // \todo IMPLEMENT Need a field selection mechanism and a field naming
00923   // (ioss_name -> stk::name)  For now, select all and give the
00924   // stk field the same name as the ioss field.
00925 
00926   // Skip the attribute field that is named "attribute"
00927   if (*I == "attribute" && names.size() > 1)
00928     continue;
00929 
00930   // \todo IMPLEMENT Need to determine whether these are
00931   // multi-state fields or constant, or interpolated, or ...
00932   Ioss::Field io_field = entity->get_field(*I);
00933   declare_ioss_field(meta, part_type, part, io_field, use_cartesian_for_scalar);
00934   }
00935 }
00936 
00937 template <typename INT>
00938 void get_entity_list(Ioss::GroupingEntity *io_entity,
00939                      stk::mesh::EntityRank part_type,
00940                      const stk::mesh::BulkData &bulk,
00941                      std::vector<stk::mesh::Entity*> &entities, INT /*dummy*/)
00942 {
00943   std::vector<INT> ids ;
00944   io_entity->get_field_data("ids", ids);
00945   
00946   size_t count = ids.size();
00947   entities.reserve(count);
00948 
00949   for(size_t i=0; i<count; ++i) {
00950     entities.push_back(bulk.get_entity( part_type, ids[i] ));
00951   }
00952 }
00953 
00954 void get_entity_list(Ioss::GroupingEntity *io_entity,
00955                      stk::mesh::EntityRank part_type,
00956                      const stk::mesh::BulkData &bulk,
00957                      std::vector<stk::mesh::Entity*> &entities)
00958 {
00959   if (db_api_int_size(io_entity) == 4) {
00960     get_entity_list(io_entity, part_type, bulk, entities, (int)0);
00961   } else {
00962     get_entity_list(io_entity, part_type, bulk, entities, (int64_t)0);
00963   }
00964 }
00965 
00966 void field_data_from_ioss(const stk::mesh::FieldBase *field,
00967                           std::vector<stk::mesh::Entity*> &entities,
00968                           Ioss::GroupingEntity *io_entity,
00969                           const std::string &io_fld_name)
00970 {
00973 
00974   if (field != NULL && io_entity->field_exists(io_fld_name)) {
00975   const Ioss::Field &io_field = io_entity->get_fieldref(io_fld_name);
00976   if (field->type_is<double>()) {
00977     internal_field_data_from_ioss(io_field, field, entities, io_entity,
00978                                     static_cast<double>(1.0));
00979   } else if (field->type_is<int>()) {
00980     // Make sure the IO field type matches the STK field type.
00981     // By default, all IO fields are created of type 'double'
00982     if (db_api_int_size(io_entity) == 4) {
00983       if (io_field.get_type() != Ioss::Field::INTEGER) {
00984         Ioss::Field &tmp = const_cast<Ioss::Field&>(io_field);
00985         tmp.reset_type(Ioss::Field::INTEGER);
00986       }
00987       internal_field_data_from_ioss(io_field, field, entities, io_entity,
00988                                       static_cast<int>(1));
00989     } else {
00990       if (io_field.get_type() != Ioss::Field::INT64) {
00991         Ioss::Field &tmp = const_cast<Ioss::Field&>(io_field);
00992         tmp.reset_type(Ioss::Field::INT64);
00993       }
00994       internal_field_data_from_ioss(io_field, field, entities, io_entity,
00995                                       static_cast<int64_t>(1));
00996     }
00997   }
00998   }
00999 }
01000 
01001 void field_data_to_ioss(const stk::mesh::FieldBase *field,
01002                         std::vector<stk::mesh::Entity*> &entities,
01003                         Ioss::GroupingEntity *io_entity,
01004                         const std::string &io_fld_name,
01005                         Ioss::Field::RoleType filter_role)
01006 {
01009 
01010   if (field != NULL && io_entity->field_exists(io_fld_name)) {
01011   const Ioss::Field &io_field = io_entity->get_fieldref(io_fld_name);
01012   if (io_field.get_role() == filter_role) {
01013     if (field->type_is<double>()) {
01014       internal_field_data_to_ioss(io_field, field, entities, io_entity,
01015                                     static_cast<double>(1.0));
01016     } else if (field->type_is<int>()) {
01017       if (io_field.get_type() != Ioss::Field::INTEGER) {
01018         Ioss::Field &tmp = const_cast<Ioss::Field&>(io_field);
01019         tmp.reset_type(Ioss::Field::INTEGER);
01020       }
01021       // FIX 64?
01022       internal_field_data_to_ioss(io_field, field, entities, io_entity,
01023                                     static_cast<int>(1));
01024     }
01025   }
01026   }
01027 }
01028 
01029 //----------------------------------------------------------------------
01030 
01031 // ----------------------------------------------------------------------
01032 // Returns true if 'entity' should be a 'part' in the analysis mesh.
01033 // Returns false if the application is only using a subset of the
01034 // database entities and this entity is not to be used. The
01035 // "omitted" property is set by the application during parsing or
01036 // pre-mesh reading time.
01037 bool include_entity(const Ioss::GroupingEntity *entity)
01038 {
01039   assert(entity);
01040 
01041   // Check whether entity has "omitted" property...
01042   bool omitted = (entity->property_exists("omitted")) &&
01043   (entity->get_property("omitted").get_int() == 1);
01044 
01045   return !omitted;
01046 }
01047 
01048 namespace {
01049 
01050 void define_side_block(stk::mesh::Part &part,
01051                        Ioss::SideSet *sset,
01052                        stk::mesh::EntityRank type,
01053                        size_t side_count, int spatial_dimension)
01054 {
01055   const stk::mesh::EntityRank siderank = side_rank(mesh::MetaData::get(part));
01056   const stk::mesh::EntityRank edgerank = edge_rank(mesh::MetaData::get(part));
01057   ThrowRequire(type == siderank || type == edgerank);
01058 
01059   const CellTopologyData *const side_topology =   stk::io::get_cell_topology(part) ?
01060     stk::io::get_cell_topology(part) :
01061     stk::mesh::fem::FEMMetaData::get(part).get_cell_topology(part).getCellTopologyData();
01062 
01063   if (side_topology == NULL) {
01064     std::ostringstream msg ;
01065     msg << " INTERNAL_ERROR: Part " << part.name() << " returned NULL from get_cell_topology()";
01066     throw std::runtime_error( msg.str() );
01067   }
01068 
01069   std::string io_topo = map_topology_cell_to_ioss(side_topology, spatial_dimension);
01070   std::string element_topo_name = "unknown";
01071 
01072   // Get sideblock parent element topology quantities...
01073   // Try to decode from part name...
01074   std::vector<std::string> tokens;
01075   stk::util::tokenize(part.name(), "_", tokens);
01076   if (tokens.size() >= 4) {
01077     // Name of form: "name_eltopo_sidetopo_id" or
01078     //               "name_block_id_sidetopo_id"
01079     // "name" is typically "surface".
01080     const Ioss::ElementTopology *element_topo = Ioss::ElementTopology::factory(tokens[tokens.size()-3], true);
01081     if (element_topo != NULL) {
01082       element_topo_name = element_topo->name();
01083     }
01084   }
01085 
01086   Ioss::SideBlock *side_block = new Ioss::SideBlock( sset->get_database() ,
01087                                                      part.name() ,
01088                                                      io_topo, element_topo_name, side_count);
01089   assert(sset->get_side_block(part.name()) == NULL);
01090   sset->add(side_block);
01091 
01092   const mesh::Field<double, mesh::ElementNode> *df = get_distribution_factor_field(part);
01093   if (df != NULL) {
01094     int nodes_per_side = side_topology->node_count;
01095     std::string storage_type = "Real[";
01096     storage_type += Ioss::Utils::to_string(nodes_per_side);
01097     storage_type += "]";
01098     side_block->field_add(Ioss::Field("distribution_factors", Ioss::Field::REAL, storage_type,
01099                                       Ioss::Field::MESH, side_count));
01100   }
01101 
01102   // Add the attribute fields.
01103   ioss_add_fields(part, part_primary_entity_rank(part), side_block, Ioss::Field::ATTRIBUTE);
01104 }
01105 
01106 void define_side_blocks(stk::mesh::Part &part,
01107                         const stk::mesh::BulkData &bulk_data,
01108                         Ioss::SideSet *sset,
01109                         stk::mesh::EntityRank type,
01110                         int spatial_dimension,
01111                         const stk::mesh::Selector *anded_selector)
01112 {
01113   mesh::MetaData & meta = mesh::MetaData::get(part);
01114   ThrowRequire(type == face_rank(meta) || type == edge_rank(meta));
01115 
01116   const stk::mesh::PartVector &blocks = part.subsets();
01117   if (blocks.size() > 0) {
01118     for (size_t j = 0; j < blocks.size(); j++) {
01119       mesh::Part & side_block_part = *blocks[j];
01120       stk::mesh::EntityRank side_rank = side_block_part.primary_entity_rank();
01121       mesh::Selector selector = meta.locally_owned_part() & side_block_part;
01122       if (anded_selector) selector &= *anded_selector;
01123 
01124       size_t num_side = count_selected_entities(selector, bulk_data.buckets(side_rank));
01125 
01126       define_side_block(side_block_part, sset, side_rank, num_side, spatial_dimension);
01127     }
01128   } else {
01129     stk::mesh::EntityRank side_rank = part.primary_entity_rank();
01130     mesh::Selector selector = meta.locally_owned_part() & part;
01131     if (anded_selector) selector &= *anded_selector;
01132     size_t num_side = count_selected_entities(selector, bulk_data.buckets(side_rank));
01133     define_side_block(part, sset, side_rank, num_side, spatial_dimension);
01134   }
01135 }
01136 
01137 //----------------------------------------------------------------------
01138 //----------------------------------------------------------------------
01139 //----------------------------------------------------------------------
01140 
01141 void define_node_block(stk::mesh::Part &part,
01142                        const stk::mesh::BulkData &bulk,
01143                        Ioss::Region &io_region,
01144                        const stk::mesh::Selector *anded_selector)
01145 {
01146   //--------------------------------
01147   // Set the spatial dimension:
01148   mesh::MetaData & meta = mesh::MetaData::get(part);
01149 
01154   mesh::Field<double, mesh::Cartesian> *coord_field =
01155     meta.get_field<stk::mesh::Field<double, mesh::Cartesian> >(std::string("coordinates"));
01156   assert(coord_field != NULL);
01157   const mesh::FieldBase::Restriction &res = coord_field->restriction(node_rank(meta), part);
01158 
01162   const int spatial_dim = res.dimension() ;
01163   io_region.property_add( Ioss::Property("spatial_dimension", spatial_dim));
01164 
01165   //--------------------------------
01166   // Create the special universal node block:
01167 
01168   mesh::Selector selector = meta.locally_owned_part() | meta.globally_shared_part();
01169   if (anded_selector) selector &= *anded_selector;
01170 
01171   size_t num_nodes = count_selected_entities(selector, bulk.buckets(node_rank(meta)));
01172 
01173   const std::string name("nodeblock_1");
01174 
01175   Ioss::NodeBlock * const nb = new Ioss::NodeBlock(io_region.get_database(),
01176                                                    name, num_nodes, spatial_dim);
01177   io_region.add( nb );
01178 
01179   // Add the attribute fields.
01180   ioss_add_fields(part, part_primary_entity_rank(part), nb, Ioss::Field::ATTRIBUTE);
01181 }
01182 
01183 
01184 void define_element_block(stk::mesh::Part &part,
01185                           const stk::mesh::BulkData &bulk,
01186                           Ioss::Region &io_region,
01187                           const stk::mesh::Selector *anded_selector)
01188 {
01189 
01190   mesh::MetaData & meta = mesh::MetaData::get(part);
01191   const stk::mesh::EntityRank elem_rank = element_rank(meta);
01192 
01193   const CellTopologyData * const cell_top =
01194     stk::io::get_cell_topology(part) ?
01195     stk::io::get_cell_topology(part) :
01196     stk::mesh::fem::FEMMetaData::get(part).get_cell_topology(part).getCellTopologyData();
01197 
01198   if (cell_top == NULL) {
01199     std::ostringstream msg ;
01200     msg << " INTERNAL_ERROR: Part " << part.name() << " returned NULL from get_cell_topology()";
01201     throw std::runtime_error( msg.str() );
01202   }
01203 
01204   mesh::Selector selector = meta.locally_owned_part() & part;
01205   if (anded_selector) selector &= *anded_selector;
01206   const size_t num_elems = count_selected_entities( selector, bulk.buckets(elem_rank));
01207 
01208   int spatial_dim = io_region.get_property("spatial_dimension").get_int();
01209 
01210   // Defer the counting of attributes until after we define the
01211   // element block so we can count them as we add them as fields to
01212   // the element block
01213   Ioss::ElementBlock *eb = new Ioss::ElementBlock(io_region.get_database() ,
01214                                                   part.name() ,
01215                                                   map_topology_cell_to_ioss(cell_top, spatial_dim) ,
01216                                                   num_elems);
01217   io_region.add(eb);
01218 
01219   // Add the attribute fields.
01220   ioss_add_fields(part, part_primary_entity_rank(part), eb, Ioss::Field::ATTRIBUTE);
01221 }
01222 
01223 void define_side_set(stk::mesh::Part &part,
01224                      const stk::mesh::BulkData &bulk,
01225                      Ioss::Region &io_region,
01226                      const stk::mesh::Selector *anded_selector)
01227 {
01228   const stk::mesh::EntityRank si_rank = side_rank(mesh::MetaData::get(part));
01229 
01230   bool create_sideset = true;
01231   if (part.subsets().empty()) {
01232     // Only define a sideset for this part if its superset part is
01233     // not a side-containing part..  (i.e., this part is not a subset part
01234     // in a surface...)
01235     const stk::mesh::PartVector &supersets = part.supersets();
01236     for (size_t i=0; i < supersets.size(); i++) {
01237       if (is_part_io_part(*supersets[i]) && supersets[i]->primary_entity_rank() == si_rank) {
01238         create_sideset = false;
01239         break;
01240       }
01241     }
01242   }
01243   if (create_sideset) {
01244     Ioss::SideSet * const ss = new Ioss::SideSet(io_region.get_database(), part.name());
01245 
01246     io_region.add(ss);
01247     int spatial_dim = io_region.get_property("spatial_dimension").get_int();
01248     define_side_blocks(part, bulk, ss, si_rank, spatial_dim, anded_selector);
01249   }
01250 }
01251 
01252 void define_node_set(stk::mesh::Part &part,
01253                      const stk::mesh::BulkData &bulk,
01254                      Ioss::Region &io_region,
01255                      const stk::mesh::Selector *anded_selector)
01256 {
01257   mesh::MetaData & meta = mesh::MetaData::get(part);
01258 
01259   mesh::Selector selector = ( meta.locally_owned_part() | meta.globally_shared_part() ) & part;
01260   if (anded_selector) selector &= *anded_selector;
01261 
01262   const size_t num_nodes =
01263     count_selected_entities(selector, bulk.buckets(node_rank(meta)));
01264 
01265   Ioss::NodeSet * const ns =
01266     new Ioss::NodeSet( io_region.get_database(), part.name(), num_nodes);
01267   io_region.add(ns);
01268 
01269   // Add the attribute fields.
01270   ioss_add_fields(part, part_primary_entity_rank(part), ns, Ioss::Field::ATTRIBUTE);
01271 }
01272 } // namespace <blank>
01273 
01274 struct part_compare {
01275   bool operator() (stk::mesh::Part *i, stk::mesh::Part *j) { return (i->name() < j->name()); }
01276 };
01277 
01278 void define_output_db(Ioss::Region & io_region ,
01279                       const mesh::BulkData &bulk_data,
01280                       const Ioss::Region *input_region,
01281                       const stk::mesh::Selector *anded_selector,
01282                       const bool sort_stk_parts)
01283 {
01284   const mesh::MetaData & meta_data = mesh::MetaData::get(bulk_data);
01285 
01286   const stk::mesh::EntityRank no_rank = node_rank(meta_data);
01287   const stk::mesh::EntityRank el_rank = element_rank(meta_data);
01288   const stk::mesh::EntityRank fa_rank = face_rank(meta_data);
01289   const stk::mesh::EntityRank ed_rank = edge_rank(meta_data);
01290 
01291   io_region.begin_mode( Ioss::STATE_DEFINE_MODEL );
01292 
01293   define_node_block(meta_data.universal_part(), bulk_data, io_region, anded_selector);
01294 
01295   // All parts of the meta data:
01296   //const mesh::PartVector & all_parts = meta_data.get_parts();
01297   const mesh::PartVector & all_parts_unsorted = meta_data.get_parts();
01298 
01299   // sort parts so they go out the same on all processors (srk: this was induced by streaming refine)
01300   mesh::PartVector all_parts = all_parts_unsorted;
01301   if (sort_stk_parts)
01302     std::sort(all_parts.begin(), all_parts.end(), part_compare());
01303 
01304   for (mesh::PartVector::const_iterator i = all_parts.begin();
01305      i != all_parts.end(); ++i) {
01306 
01307   mesh::Part * const part = *i ;
01308 
01309   if (is_part_io_part(*part)) {
01310     if (invalid_rank(part->primary_entity_rank()))
01311       continue;
01312       else if (part->primary_entity_rank() == no_rank)
01313       define_node_set(*part, bulk_data, io_region, anded_selector);
01314       else if (part->primary_entity_rank() == el_rank)
01315       define_element_block(*part, bulk_data, io_region, anded_selector);
01316       else if (part->primary_entity_rank() == fa_rank)
01317       define_side_set(*part, bulk_data, io_region, anded_selector);
01318       else if (part->primary_entity_rank() == ed_rank)
01319       define_side_set(*part, bulk_data, io_region, anded_selector);
01320   }
01321   }
01322 
01323   if (input_region != NULL)
01324   io_region.synchronize_id_and_name(input_region, true);
01325 
01326   // for streaming refinement, each "pseudo-processor" doesn't know about others, so we pick a sort order
01327   //   and use it for all pseudo-procs - the original_block_order property is used to set the order
01328   //   on all procs.
01329   if (sort_stk_parts)
01330   {
01331     int offset=0;
01332     for (mesh::PartVector::const_iterator i = all_parts.begin();
01333          i != all_parts.end(); ++i) {
01334 
01335       mesh::Part * const part = *i ;
01336 
01337       if (is_part_io_part(*part)) {
01338         if (invalid_rank(part->primary_entity_rank()))
01339           continue;
01340         else if (part->primary_entity_rank() == el_rank)
01341           {
01342             Ioss::GroupingEntity *element_block = io_region.get_entity(part->name());
01343             if (element_block)
01344               {
01345                 if (element_block->property_exists("original_block_order")) {
01346                   element_block->property_erase("original_block_order");
01347                 }
01348                 element_block->property_add(Ioss::Property("original_block_order", offset));
01349                 ++offset;
01350               }
01351           }
01352       }
01353     }
01354   }
01355 
01356   io_region.end_mode( Ioss::STATE_DEFINE_MODEL );
01357 }
01358 
01359 //----------------------------------------------------------------------
01360 
01361 namespace {
01362 
01363 size_t get_entities(stk::mesh::Part &part,
01364                     const stk::mesh::BulkData &bulk,
01365                     std::vector<mesh::Entity*> &entities,
01366                     bool include_shared,
01367                     const stk::mesh::Selector *anded_selector)
01368 {
01369   mesh::MetaData & meta = mesh::MetaData::get(part);
01370   mesh::EntityRank type = part_primary_entity_rank(part);
01371   if (invalid_rank(type))
01372     type = node_rank(meta);
01373 
01374   mesh::Selector own_share = meta.locally_owned_part();
01375   if (include_shared)
01376     own_share |= meta.globally_shared_part();
01377 
01378   mesh::Selector selector = part & own_share;
01379   if (anded_selector) selector &= *anded_selector;
01380 
01381   get_selected_entities(selector, bulk.buckets(type), entities);
01382   return entities.size();
01383 }
01384 
01385 
01386 template <typename INT>
01387 void write_side_data_to_ioss( Ioss::GroupingEntity & io ,
01388                               mesh::Part * const part ,
01389                               const mesh::BulkData & bulk_data,
01390                               const stk::mesh::Selector *anded_selector, INT /*dummy*/ )
01391 {
01392   //std::cout << "tmp write_side_data_to_ioss part= " << part->name() << std::endl;
01393   const mesh::MetaData & meta_data = mesh::MetaData::get(*part);
01394 
01395   std::vector<mesh::Entity *> sides ;
01396   size_t num_sides = get_entities(*part, bulk_data, sides, false, anded_selector);
01397 
01398   std::vector<INT> elem_side_ids; elem_side_ids.reserve(num_sides*2);
01399 
01400   stk::mesh::EntityRank elem_rank = element_rank(meta_data);
01401   for(size_t i=0; i<num_sides; ++i) {
01402 
01403     const mesh::Entity &side = *sides[i] ;
01404     const mesh::PairIterRelation side_elem = side.relations( elem_rank );
01405 
01406     // Which element to use?
01407     // Any locally owned element that has the "correct" orientation
01408 
01409     const size_t num_side_elem = side_elem.size();
01410 
01411     const mesh::Relation *rel = NULL ;
01412 
01413     for ( size_t j = 0 ; j < num_side_elem && ! rel ; ++j ) {
01414       const mesh::Entity & elem = *side_elem[j].entity();
01415 
01416       if ( elem.bucket().member( meta_data.locally_owned_part() ) &&
01417            (num_side_elem == 1 || stk::mesh::fem::element_side_polarity(elem, side, side_elem[j].identifier())) ) {
01418         rel = &side_elem[j];
01419       }
01420     }
01421 
01422     if (rel == NULL) { // no suitable element found
01423       std::ostringstream oss;
01424       oss << "ERROR, no suitable element found";
01425       throw std::runtime_error(oss.str());
01426     }
01427 
01428     elem_side_ids.push_back(rel->entity()->identifier());
01429     elem_side_ids.push_back(rel->identifier() + 1) ; // Ioss is 1-based, mesh is 0-based.
01430   }
01431 
01432   const size_t num_side_written = io.put_field_data("element_side",elem_side_ids);
01433 
01434   if ( num_sides != num_side_written ) {
01435     std::ostringstream msg ;
01436 
01437     msg << "stk::io::write_side_data_to_ioss FAILED for " ;
01438     msg << io.name();
01439     msg << " in Ioss::GroupingEntity::put_field_data:" ;
01440     msg << " num_sides = " << num_sides ;
01441     msg << " , num_side_written = " << num_side_written ;
01442     throw std::runtime_error( msg.str() );
01443   }
01444 
01445   const mesh::Field<double, mesh::ElementNode> *df = get_distribution_factor_field(*part);
01446   if (df != NULL) {
01447     field_data_to_ioss(df, sides, &io, "distribution_factors", Ioss::Field::MESH);
01448   }
01449 
01450   const std::vector<mesh::FieldBase *> &fields = meta_data.get_fields();
01451   std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
01452   while (I != fields.end()) {
01453     const mesh::FieldBase *f = *I ; ++I ;
01454     const Ioss::Field::RoleType *role = stk::io::get_field_role(*f);
01455     if (role != NULL && *role == Ioss::Field::ATTRIBUTE) {
01456       stk::io::field_data_to_ioss(f, sides, &io, f->name(), Ioss::Field::ATTRIBUTE);
01457     }
01458   }
01459 }
01460 
01461 //----------------------------------------------------------------------
01462 template <typename INT>
01463 void output_node_block(Ioss::NodeBlock &nb,
01464                        stk::mesh::Part &part,
01465                        const stk::mesh::BulkData &bulk,
01466                        const stk::mesh::Selector *anded_selector, INT /*dummy*/)
01467 {
01468   //----------------------------------
01469   // Exactly one node block to obtain the nodal coordinates and ids:
01470   // Note that the "ids" field of the nodes needs to be written
01471   // before any other bulk data that uses node ids since it sets up
01472   // the global->local mapping of nodes for the output database.
01473   // Similarly for the element "ids" field related to bulk data
01474   // using element ids.
01475   std::vector<mesh::Entity *> nodes ;
01476   size_t num_nodes = get_entities(part, bulk, nodes, true, anded_selector);
01477 
01478   std::vector<INT> node_ids; node_ids.reserve(num_nodes);
01479   for(size_t i=0; i<num_nodes; ++i) {
01480     const mesh::Entity & node = * nodes[i] ;
01481     node_ids.push_back(node.identifier());
01482   }
01483 
01484   size_t num_ids_written = nb.put_field_data("ids", node_ids);
01485   if ( num_nodes != num_ids_written) {
01486     std::ostringstream msg ;
01487     msg << " FAILED in Ioss::NodeBlock::put_field_data:" ;
01488     msg << " num_nodes = " << num_nodes ;
01489     msg << " , num_ids_written = " << num_ids_written ;
01490     throw std::runtime_error( msg.str() );
01491   }
01492 
01497   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01498   mesh::Field<double, mesh::Cartesian> *coord_field =
01499     meta_data.get_field<stk::mesh::Field<double, mesh::Cartesian> >(std::string("coordinates"));
01500   assert(coord_field != NULL);
01501   field_data_to_ioss(coord_field, nodes, &nb, "mesh_model_coordinates", Ioss::Field::MESH);
01502 
01503   const std::vector<mesh::FieldBase *> &fields = meta_data.get_fields();
01504   std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
01505   while (I != fields.end()) {
01506     const mesh::FieldBase *f = *I ; ++I ;
01507     if (stk::io::is_valid_part_field(f, part_primary_entity_rank(part), part,
01508              meta_data.universal_part(), Ioss::Field::ATTRIBUTE, false)) {
01509       stk::io::field_data_to_ioss(f, nodes, &nb, f->name(), Ioss::Field::ATTRIBUTE);
01510     }
01511   }
01512 }
01513 
01514 template <typename INT>
01515 void output_element_block(Ioss::ElementBlock *block,
01516                           const stk::mesh::BulkData &bulk,
01517                           const stk::mesh::Selector *anded_selector, INT /*dummy*/)
01518 {
01519   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01520   const std::string& name = block->name();
01521   mesh::Part* part = meta_data.get_part(name);
01522 
01523   assert(part != NULL);
01524   std::vector<mesh::Entity *> elements;
01525   size_t num_elems = get_entities(*part, bulk, elements, false, anded_selector);
01526 
01527   const CellTopologyData * cell_topo =
01528     stk::io::get_cell_topology(*part) ?
01529     stk::io::get_cell_topology(*part) :
01530     stk::mesh::fem::FEMMetaData::get(*part).get_cell_topology(*part).getCellTopologyData();
01531   if (cell_topo == NULL) {
01532     std::ostringstream msg ;
01533     msg << " INTERNAL_ERROR: Part " << part->name() << " returned NULL from get_cell_topology()";
01534     throw std::runtime_error( msg.str() );
01535   }
01536   size_t nodes_per_elem = cell_topo->node_count;
01537 
01538   std::vector<INT> elem_ids; elem_ids.reserve(num_elems);
01539   std::vector<INT> connectivity; connectivity.reserve(num_elems*nodes_per_elem);
01540 
01541   stk::mesh::EntityRank no_rank = node_rank(meta_data);
01542   for (size_t i = 0; i < num_elems; ++i) {
01543 
01544     elem_ids.push_back(elements[i]->identifier());
01545 
01546     const mesh::PairIterRelation elem_nodes = elements[i]->relations(no_rank);
01547 
01548     for (size_t j = 0; j < nodes_per_elem; ++j) {
01549       connectivity.push_back(elem_nodes[j].entity()->identifier());
01550     }
01551   }
01552 
01553   const size_t num_ids_written = block->put_field_data("ids", elem_ids);
01554   const size_t num_con_written = block->put_field_data("connectivity", connectivity);
01555 
01556   if ( num_elems != num_ids_written || num_elems != num_con_written ) {
01557     std::ostringstream msg ;
01558     msg << " FAILED in Ioss::ElementBlock::put_field_data:" << std::endl ;
01559     msg << "  num_elems = " << num_elems << std::endl ;
01560     msg << "  num_ids_written = " << num_ids_written << std::endl ;
01561     msg << "  num_connectivity_written = " << num_con_written << std::endl ;
01562     throw std::runtime_error( msg.str() );
01563   }
01564 
01565   stk::mesh::EntityRank elem_rank = element_rank(meta_data);
01566   const std::vector<mesh::FieldBase *> &fields = meta_data.get_fields();
01567   std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
01568   while (I != fields.end()) {
01569     const mesh::FieldBase *f = *I ; ++I ;
01570     const Ioss::Field::RoleType *role = stk::io::get_field_role(*f);
01571     if (role != NULL && *role == Ioss::Field::ATTRIBUTE) {
01572       const mesh::FieldBase::Restriction &res = f->restriction(elem_rank, *part);
01573       if (res.dimension() > 0) {
01574         stk::io::field_data_to_ioss(f, elements, block, f->name(), Ioss::Field::ATTRIBUTE);
01575       }
01576     }
01577   }
01578 }
01579 
01580 template <typename INT>
01581 void output_node_set(Ioss::NodeSet *ns, const stk::mesh::BulkData &bulk,
01582                      const stk::mesh::Selector *anded_selector, INT /*dummy*/)
01583 {
01584   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01585   const std::string& name = ns->name();
01586   stk::mesh::Part* part = meta_data.get_part(name);
01587   assert(part != NULL);
01588 
01589   std::vector<stk::mesh::Entity *> nodes ;
01590   size_t num_nodes = get_entities(*part, bulk, nodes, true, anded_selector);
01591 
01592   std::vector<INT> node_ids; node_ids.reserve(num_nodes);
01593   for(size_t i=0; i<num_nodes; ++i) {
01594     const stk::mesh::Entity & node = * nodes[i] ;
01595     node_ids.push_back(node.identifier());
01596   }
01597 
01598   size_t num_ids_written = ns->put_field_data("ids", node_ids);
01599   if ( num_nodes != num_ids_written ) {
01600     std::ostringstream msg ;
01601     msg << " FAILED in Ioss::NodeSet::put_field_data:"
01602         << " num_nodes = " << num_nodes
01603         << ", num_ids_written = " << num_ids_written;
01604     throw std::runtime_error( msg.str() );
01605   }
01606 
01607   stk::mesh::Field<double> *df_field =
01608     meta_data.get_field<stk::mesh::Field<double> >("distribution_factors");
01609   if (df_field != NULL) {
01610     stk::io::field_data_to_ioss(df_field, nodes, ns, "distribution_factors", Ioss::Field::MESH);
01611   }
01612 
01613   const std::vector<mesh::FieldBase *> &fields = meta_data.get_fields();
01614   std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
01615   while (I != fields.end()) {
01616     const mesh::FieldBase *f = *I ; ++I ;
01617     const Ioss::Field::RoleType *role = stk::io::get_field_role(*f);
01618     if (role != NULL && *role == Ioss::Field::ATTRIBUTE) {
01619       const mesh::FieldBase::Restriction &res = f->restriction(0, *part);
01620       if (res.dimension() > 0) {
01621         stk::io::field_data_to_ioss(f, nodes, ns, f->name(), Ioss::Field::ATTRIBUTE);
01622       }
01623     }
01624   }
01625 }
01626 
01627 template <typename INT>
01628 void output_side_set(Ioss::SideSet *ss,
01629                      const stk::mesh::BulkData &bulk,
01630                      const stk::mesh::Selector *anded_selector, INT dummy)
01631 {
01632   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01633   size_t block_count = ss->block_count();
01634   for (size_t i=0; i < block_count; i++) {
01635     Ioss::SideBlock *block = ss->get_block(i);
01636     if (stk::io::include_entity(block)) {
01637       stk::mesh::Part * const part = meta_data.get_part(block->name());
01638       stk::io::write_side_data_to_ioss(*block, part, bulk, anded_selector, dummy);
01639     }
01640   }
01641 }
01642 
01643 } // namespace <blank>
01644 
01645 void write_output_db(Ioss::Region& io_region,
01646                      const stk::mesh::BulkData& bulk,
01647                      const stk::mesh::Selector *anded_selector)
01648 {
01649   const stk::mesh::MetaData & meta = mesh::MetaData::get(bulk);
01650 
01651   bool ints64bit = db_api_int_size(&io_region) == 8;
01652 
01653   io_region.begin_mode( Ioss::STATE_MODEL );
01654 
01655   int64_t z64 = 0;
01656   int     z32 = 0;
01657 
01658   Ioss::NodeBlock & nb = *io_region.get_node_blocks()[0];
01659   if (ints64bit)
01660   output_node_block(nb, meta.universal_part(), bulk, anded_selector, z64);
01661   else
01662   output_node_block(nb, meta.universal_part(), bulk, anded_selector, z32);
01663 
01664   //----------------------------------
01665   const Ioss::ElementBlockContainer& elem_blocks = io_region.get_element_blocks();
01666   for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
01667     it != elem_blocks.end(); ++it) {
01668   if (ints64bit)
01669     output_element_block(*it, bulk, anded_selector, z64);
01670   else
01671     output_element_block(*it, bulk, anded_selector, z32);
01672   }
01673 
01674   //----------------------------------
01675   const Ioss::NodeSetContainer& node_sets = io_region.get_nodesets();
01676   for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
01677     it != node_sets.end(); ++it) {
01678   if (ints64bit)
01679     output_node_set(*it, bulk, anded_selector, z64);
01680   else
01681     output_node_set(*it, bulk, anded_selector, z32);
01682   }
01683 
01684   //----------------------------------
01685   const Ioss::SideSetContainer& side_sets = io_region.get_sidesets();
01686   for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
01687     it != side_sets.end(); ++it) {
01688   if (ints64bit)
01689     output_side_set(*it, bulk, anded_selector, z64);
01690   else
01691     output_side_set(*it, bulk, anded_selector, z32);
01692   }
01693 
01694   io_region.end_mode( Ioss::STATE_MODEL );
01695 }
01696 
01697 //----------------------------------------------------------------------
01698 bool is_part_io_part(stk::mesh::Part &part)
01699 {
01700   return NULL != part.attribute<Ioss::GroupingEntity>();
01701 }
01702 
01703 const stk::mesh::Field<double, stk::mesh::ElementNode> *get_distribution_factor_field(const stk::mesh::Part &p)
01704 {
01705   return p.attribute<stk::mesh::Field<double, stk::mesh::ElementNode> >();
01706 }
01707 
01708 void set_distribution_factor_field(stk::mesh::Part &p,
01709                                    const stk::mesh::Field<double, stk::mesh::ElementNode> &df_field)
01710 {
01711   stk::mesh::MetaData &m = mesh::MetaData::get(p);
01712   m.declare_attribute_no_delete(p,&df_field);
01713 }
01714 
01715 const Ioss::Field::RoleType* get_field_role(const stk::mesh::FieldBase &f)
01716 {
01717   return f.attribute<Ioss::Field::RoleType>();
01718 }
01719 
01720 void set_field_role(stk::mesh::FieldBase &f, const Ioss::Field::RoleType &role)
01721 {
01722   Ioss::Field::RoleType *my_role = new Ioss::Field::RoleType(role);
01723   stk::mesh::MetaData &m = mesh::MetaData::get(f);
01724   const Ioss::Field::RoleType *check = m.declare_attribute_with_delete(f, my_role);
01725   if ( check != my_role ) {
01726     if (*check != *my_role) {
01727       std::ostringstream msg ;
01728     msg << " FAILED in IossBridge -- set_field_role:"
01729         << " The role type had already been set to " << *check
01730         << ", so it is not possible to change it to " << *my_role;
01731     throw std::runtime_error( msg.str() );
01732     }
01733     delete my_role;
01734   }
01735 }
01736 
01737 }//namespace io
01738 }//namespace stk
01739 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines