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