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(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   Ioss::SideSet *sset = dynamic_cast<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   Ioss::SideBlock *sblk = dynamic_cast<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     int 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     int 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(int 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     int 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(int j=0; j<field_component_count; ++j) {
00378     io_field_data[i*field_component_count+j] = fld_data[j];
00379   }
00380       } else {
00381   for(int 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     bool invalid_rank(stk::mesh::EntityRank rank)
00406     {
00407       return rank == mesh::InvalidEntityRank;
00408     }
00409 
00410     stk::mesh::EntityRank part_primary_entity_rank(const stk::mesh::Part &part)
00411     {
00412       if (mesh::MetaData::get(part).universal_part() == part) {
00413         return stk::mesh::fem::FEMMetaData::NODE_RANK;
00414       }
00415       else {
00416         return part.primary_entity_rank();
00417       }
00418     }
00419 
00420     stk::mesh::EntityRank element_rank(const stk::mesh::MetaData &meta)
00421     {
00422       return stk::mesh::fem::FEMMetaData::get(meta).element_rank();
00423     }
00424 
00425     stk::mesh::EntityRank side_rank(const stk::mesh::MetaData &meta)
00426     {
00427       return stk::mesh::fem::FEMMetaData::get(meta).side_rank();
00428     }
00429 
00430     stk::mesh::EntityRank face_rank(const stk::mesh::MetaData &meta)
00431     {
00432       return stk::mesh::fem::FEMMetaData::get(meta).face_rank();
00433     }
00434 
00435     stk::mesh::EntityRank edge_rank(const stk::mesh::MetaData &meta)
00436     {
00437       return stk::mesh::fem::FEMMetaData::get(meta).edge_rank();
00438     }
00439 
00440     stk::mesh::EntityRank node_rank(const stk::mesh::MetaData& meta)
00441     {
00442       return stk::mesh::fem::FEMMetaData::NODE_RANK;
00443     }
00444 
00445     void set_cell_topology(stk::mesh::Part &part, const CellTopologyData * const cell_topology)
00446     {
00447       stk::mesh::fem::set_cell_topology(part, cell_topology);
00448     }
00449 
00450     const CellTopologyData *get_cell_topology(const stk::mesh::Part &part)
00451     {
00452       return stk::mesh::fem::FEMMetaData::get(part).get_cell_topology(part).getCellTopologyData();
00453     }
00454 
00455     void initialize_spatial_dimension(stk::mesh::fem::FEMMetaData & fem_meta, size_t spatial_dimension,
00456               const std::vector<std::string> &entity_rank_names)
00457     {
00458       if (!fem_meta.is_FEM_initialized() ) {
00459         fem_meta.FEM_initialize(spatial_dimension, entity_rank_names);
00460       }
00461     }
00462 
00463     void initialize_spatial_dimension(stk::mesh::MetaData &meta, size_t spatial_dimension,
00464               const std::vector<std::string> &entity_rank_names)
00465     {
00466       stk::mesh::fem::FEMMetaData & fem_meta = stk::mesh::fem::FEMMetaData::get(meta);
00467       initialize_spatial_dimension(fem_meta, spatial_dimension, entity_rank_names);
00468     }
00469 
00470     void get_io_field_type(const stk::mesh::FieldBase *field,
00471          const stk::mesh::FieldRestriction &res,
00472          std::pair<std::string, Ioss::Field::BasicType> *result)
00473     {
00474       static const std::string invalid("invalid");
00475       static const std::string scalar("scalar");
00476       static const std::string vector_2d("vector_2d");
00477       static const std::string vector_3d("vector_3d");
00478       static const std::string quaternion_2d("quaternion_2d");
00479       static const std::string quaternion_3d("quaternion_3d");
00480       static const std::string full_tensor_36("full_tensor_36");
00481       static const std::string full_tensor_32("full_tensor_32");
00482       static const std::string full_tensor_22("full_tensor_22");
00483       static const std::string full_tensor_16("full_tensor_16");
00484       static const std::string full_tensor_12("full_tensor_12");
00485       static const std::string sym_tensor_33("sym_tensor_33");
00486       static const std::string sym_tensor_31("sym_tensor_31");
00487       static const std::string sym_tensor_21("sym_tensor_21");
00488       static const std::string sym_tensor_13("sym_tensor_13");
00489       static const std::string sym_tensor_11("sym_tensor_11");
00490       static const std::string sym_tensor_10("sym_tensor_10");
00491       static const std::string asym_tensor_03("asym_tensor_03");
00492       static const std::string asym_tensor_02("asym_tensor_02");
00493       static const std::string asym_tensor_01("asym_tensor_01");
00494       static const std::string matrix_22("matrix_22");
00495       static const std::string matrix_33("matrix_33");
00496 
00497       const unsigned rank = field->rank();
00498       const shards::ArrayDimTag * const * const tags = field->dimension_tags();
00499 
00500       result->second = Ioss::Field::INVALID;
00501 
00502       if ( field->type_is<double>() ) {
00503   result->second = Ioss::Field::REAL;
00504       }
00505       else if ( field->type_is<int>() ) {
00506   result->second = Ioss::Field::INTEGER;
00507       }
00508 
00509       if ( 0 == rank ) {
00510   result->first = scalar ;
00511       }
00512       else if ( 1 == rank ) {
00513   size_t num_comp = res.stride(0);
00514   if ( tags[0] == & stk::mesh::Cartesian::tag() && 1 == num_comp ) {
00515     result->first = scalar ;
00516   }
00517   else if ( tags[0] == & stk::mesh::Cartesian::tag() && 2 == num_comp ) {
00518     result->first = vector_2d ;
00519   }
00520   else if ( tags[0] == & stk::mesh::Cartesian::tag() && 3 == num_comp ) {
00521     result->first = vector_3d ;
00522   }
00523   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 9 == num_comp ) {
00524     result->first = full_tensor_36 ;
00525   }
00526   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 5 == num_comp ) {
00527     result->first = full_tensor_32 ;
00528   }
00529   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 4 == num_comp ) {
00530     result->first = full_tensor_22 ;
00531   }
00532   else if ( tags[0] == & stk::mesh::FullTensor::tag() && 3 == num_comp ) {
00533     result->first = full_tensor_12 ;
00534   }
00535   else if ( tags[0] == & stk::mesh::SymmetricTensor::tag() && 6 == num_comp ) {
00536     result->first = sym_tensor_33 ;
00537   }
00538   else if ( tags[0] == & stk::mesh::SymmetricTensor::tag() && 4 == num_comp ) {
00539     result->first = sym_tensor_31 ;
00540   }
00541   else if ( tags[0] == & stk::mesh::SymmetricTensor::tag() && 3 == num_comp ) {
00542     result->first = sym_tensor_21 ;
00543   }
00544       }
00545 
00546       if ( result->first.empty() ) {
00547   size_t num_comp = res.stride(rank-1);
00548   std::ostringstream tmp ;
00549   tmp << "Real[" << num_comp << "]" ;
00550   result->first = tmp.str();
00551       }
00552     }
00553 
00554     //----------------------------------------------------------------------
00555     const Ioss::GroupingEntity *get_associated_ioss_entity(const mesh::Part &part)
00556     {
00557       const Ioss::GroupingEntity *entity = part.attribute<Ioss::GroupingEntity>();
00558       if (!entity || entity->type() == Ioss::INVALID_TYPE) {
00559   return NULL;
00560       } else {
00561   return entity;
00562       }
00563     }
00564 
00565     void put_io_part_attribute(mesh::Part & part, Ioss::GroupingEntity *entity)
00566     {
00567       if (part.attribute<Ioss::GroupingEntity>() != NULL) {
00568   std::string msg = "stk::io::put_io_part_attribute( ";
00569   msg += part.name();
00570   msg += " ) FAILED:";
00571   msg += " io_part_attribute is already defined";
00572   throw std::runtime_error( msg );
00573       }
00574 
00575       mesh::MetaData & meta = mesh::MetaData::get(part);
00576       if (entity) {
00577   meta.declare_attribute_no_delete(part, entity);
00578       } else {
00579   Ioss::GroupingEntity *attr = new Ioss::NullEntity();
00580   meta.declare_attribute_with_delete(part, attr);
00581       }
00582     }
00583 
00584     void remove_io_part_attribute(mesh::Part & part)
00585     {
00586       const Ioss::GroupingEntity *entity = part.attribute<Ioss::GroupingEntity>();
00587       if (entity == NULL) {
00588   std::string msg = "stk::io::remove_io_part_attribute( ";
00589   msg += part.name();
00590   msg += " ) FAILED:";
00591   msg += " io_part_attribute is not defined on this part";
00592   throw std::runtime_error( msg );
00593       } else {
00594   mesh::MetaData & meta = mesh::MetaData::get(part);
00595   bool success = meta.remove_attribute(part, entity);
00596   if (!success) {
00597     std::string msg = "stk::io::remove_io_part_attribute( ";
00598     msg += part.name();
00599     msg += " ) FAILED:";
00600     msg += " meta.remove_attribute(..) returned failure.";
00601     throw std::runtime_error( msg );
00602   }
00603 
00604   if (entity->type() == Ioss::INVALID_TYPE) {
00605     delete entity;
00606   }
00607 
00608       }
00609     }
00610 
00615     bool is_valid_part_field(const stk::mesh::FieldBase *field,
00616            stk::mesh::EntityRank part_type,
00617            stk::mesh::Part &part,
00618            stk::mesh::Part &universal,
00619            Ioss::Field::RoleType filter_role,
00620            bool add_all)
00621     {
00622       const Ioss::Field::RoleType *role = stk::io::get_field_role(*field);
00623 
00624       if (!add_all && role == NULL) {
00625   return false;
00626       }
00627 
00628       if (role != NULL && *role != filter_role)
00629   return false;
00630 
00631       const stk::mesh::FieldBase::Restriction &res = field->restriction(part_type, part);
00632       if (res.dimension() > 0) {
00633   // The field exists on the current 'part'.  Now check (for
00634   // node types only) whether the 'part' is *either* the
00635   // universal_part() *or* the field *doesn't* exist on the
00636   // universal part...
00637   // Note that for "node" type parts, the IO database has a part
00638   // (nodeblock) that corresponds to the universal_part(), so
00639   // fields defined on all nodes are output on the nodeblock and
00640   // fields defined on only a subset of the parts should be
00641   // output on that subset which maps to a nodeset.  For other
00642   // part types, there is no IO entity that corresponds to a
00643   // universal part, so fields must be output on the part they
00644   // exist on.  There may be a problem if we start using element
00645   // sets ..., but wait until we get to that point; current code
00646   // works with current entity set.
00647   if (part_type != node_rank(mesh::MetaData::get(part)) || part == universal) {
00648     return true;
00649   }
00650 
00651   const stk::mesh::FieldBase::Restriction &res_universe = field->restriction(part_type, universal);
00652   if (res_universe.dimension() <= 0) {
00653     // Field exists on current part, but not on the universal
00654     // set (and this part is not the universal part)
00655     return true;
00656   }
00657       }
00658       return false;
00659     }
00660 
00664     // ========================================================================
00665     const CellTopologyData *map_topology_ioss_to_cell(const Ioss::ElementTopology *topology)
00666     {
00673 
00674       std::string name         = topology->name();
00675       int io_nodes_per_element = topology->number_nodes();
00676 
00677       const CellTopologyData *cell_topology = map_ioss_to_topology(name, io_nodes_per_element);
00678 
00679       return cell_topology;
00680     }
00681 
00682     std::string map_topology_cell_to_ioss( const CellTopologyData *cell_top,
00683              int spatial_dimension)
00684     {
00685       std::string extype = "unknown";
00686 
00687       if (cell_top == NULL)
00688   return extype;
00689 
00690       if(strcmp(cell_top->name, "super") == 0) {
00691         std::stringstream oss;
00692         oss << "super" << cell_top->node_count;
00693         return oss.str();
00694       }
00695       else if(strncasecmp(cell_top->name, "super", 5) == 0) {
00696         return cell_top->name;
00697       }
00698 
00699       switch( cell_top->key ) {
00700       case shards::Node::key : extype.assign( "node" ); break ;
00701       case shards::Particle::key :
00702   if (spatial_dimension == 2) extype = "circle1";
00703   else                        extype = "sphere1";
00704   break ;
00705 
00706       case shards::Beam<2>::key :
00707   extype = "bar2";
00708   break ;
00709 
00710       case shards::Beam<3>::key :
00711   extype = "bar3";
00712   break ;
00713 
00714       case shards::Line<2>::key :
00715   if (spatial_dimension == 2) extype = "edge2d2";
00716   else                        extype = "edge2";
00717   break ;
00718       case shards::Line<3>::key :
00719   if (spatial_dimension == 2) extype = "edge2d3";
00720   else                        extype = "edge3";
00721   break ;
00722 
00723       case shards::ShellLine<2>::key : extype.assign( "shellline2d2" ); break ;
00724       case shards::ShellLine<3>::key : extype.assign( "shellline2d3" ); break ;
00725 
00726       case shards::Triangle<3>::key :
00727   extype = "tri3";
00728   break ;
00729       case shards::Triangle<4>::key :
00730         extype = "tri4";
00731   break ;
00732       case shards::Triangle<6>::key :
00733         extype = "tri6";
00734   break ;
00735 
00736       case shards::ShellTriangle<3>::key : extype.assign( "trishell3" ); break ;
00737       case shards::ShellTriangle<6>::key : extype.assign( "trishell6" ); break ;
00738 
00739       case shards::Quadrilateral<4>::key :
00740         extype = "quad4";
00741   break ;
00742       case shards::Quadrilateral<8>::key :
00743         extype = "quad8";
00744   break ;
00745       case shards::Quadrilateral<9>::key :
00746         extype = "quad9";
00747   break ;
00748 
00749       case shards::ShellQuadrilateral<4>::key : extype.assign( "shell4" ); break ;
00750       case shards::ShellQuadrilateral<8>::key : extype.assign( "shell8" ); break ;
00751       case shards::ShellQuadrilateral<9>::key : extype.assign( "shell9" ); break ;
00752 
00753       case shards::Tetrahedron< 4>::key : extype.assign( "tetra4" ); break ;
00754       case shards::Tetrahedron<8>::key : extype.assign( "tetra8" ); break ;
00755       case shards::Tetrahedron<10>::key : extype.assign( "tetra10" ); break ;
00756 
00757       case shards::Pyramid< 5>::key : extype.assign( "pyramid5" ); break ;
00758       case shards::Pyramid<13>::key : extype.assign( "pyramid13" ); break ;
00759       case shards::Pyramid<14>::key : extype.assign( "pyramid14" ); break ;
00760 
00761       case shards::Wedge< 6>::key : extype.assign( "wedge6" ); break ;
00762       case shards::Wedge<15>::key : extype.assign( "wedge15" ); break ;
00763       case shards::Wedge<18>::key : extype.assign( "wedge18" ); break ;
00764 
00765       case shards::Hexahedron< 8>::key : extype.assign( "hex8" ); break ;
00766       case shards::Hexahedron<20>::key : extype.assign( "hex20" ); break ;
00767       case shards::Hexahedron<27>::key : extype.assign( "hex27" ); break ;
00768 
00769       default:
00770   std::ostringstream oss;
00771   oss << "stk::io::map_topology_to_ioss( '" << cell_top->name
00772       << "' ) ERROR unmapped topology" << std::endl ;
00773   throw std::runtime_error(oss.str());
00774       }
00775 
00776       return extype ;
00777     }
00778 
00779     void internal_part_processing(Ioss::GroupingEntity *entity, stk::mesh::fem::FEMMetaData &meta)
00780     {
00781       internal_part_processing(entity, meta.get_meta_data(meta));
00782     }
00783 
00784     void internal_part_processing(Ioss::EntityBlock *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::GroupingEntity *entity, stk::mesh::MetaData &meta)
00790     {
00791       if (include_entity(entity)) {
00792         stk::mesh::fem::FEMMetaData * fem_meta = const_cast<stk::mesh::fem::FEMMetaData *>(meta.get_attribute<stk::mesh::fem::FEMMetaData>());
00793   mesh::EntityRank type = get_entity_rank(entity, meta);
00794         if (fem_meta) {
00795     stk::mesh::Part & part = fem_meta->declare_part(entity->name(), type);
00796     stk::io::put_io_part_attribute(part, entity);
00797         } else {
00798     stk::mesh::Part & part = meta.declare_part(entity->name(), type);
00799     stk::io::put_io_part_attribute(part, entity);
00800         }
00801       }
00802     }
00803 
00804     void internal_part_processing(Ioss::EntityBlock *entity, stk::mesh::MetaData &meta)
00805     {
00806       if (include_entity(entity)) {
00807   mesh::EntityRank type = get_entity_rank(entity, meta);
00808         stk::mesh::fem::FEMMetaData * fem_meta = const_cast<stk::mesh::fem::FEMMetaData *>(meta.get_attribute<stk::mesh::fem::FEMMetaData>());
00809         //const stk::mesh::fem::FEMMetaData * fem_meta = meta.get_attribute<stk::mesh::fem::FEMMetaData>();
00810         stk::mesh::Part * part = NULL;
00811         if( fem_meta )
00812           part = &fem_meta->declare_part(entity->name(), type);
00813         else
00814           part = &meta.declare_part(entity->name(), type);
00815   stk::io::put_io_part_attribute(*part, entity);
00816 
00817   const Ioss::ElementTopology *topology = entity->topology();
00818   // Check spatial dimension of the element topology here so we
00819   // can issue a more meaningful error message.  If the
00820   // dimension is bad and we continue to the following calls,
00821   // there is an exception and we get unintelligible (to the
00822   // user) error messages.  Could also do a catch...
00823 
00824   if (entity->type() == Ioss::ELEMENTBLOCK) {
00825     assert(topology != NULL);
00826     if (fem_meta && (topology->spatial_dimension() < (int)fem_meta->spatial_dimension())) {
00827       // NOTE: The comparison is '<' and not '!=' since a 2D mesh
00828       // can contain a "3d" element -- a Beam is both a 2D and
00829       // 3D element...
00830 
00831       std::ostringstream msg ;
00832       msg << "\n\nERROR: Element Block " << entity->name()
00833     << " contains " << topology->name() << " elements with spatial dimension "
00834     << topology->spatial_dimension()
00835     << "\n       which does not match the spatial dimension of the model which is "
00836     << fem_meta->spatial_dimension() << "\n\n";
00837       throw std::runtime_error( msg.str() );
00838     }
00839   }
00840 
00841   const CellTopologyData * const cell_topology = map_topology_ioss_to_cell(topology);
00846 
00847         if (cell_topology != NULL) {
00848           if( fem_meta )
00849           {
00850             const stk::mesh::fem::CellTopology cell_topo(cell_topology);
00851             stk::mesh::fem::set_cell_topology(*part, cell_topo);
00852           }
00853           stk::io::set_cell_topology(*part, cell_topology);
00854         } else {
00856   }
00857   stk::io::define_io_fields(entity, Ioss::Field::ATTRIBUTE, *part, type);
00858       }
00859     }
00860 
00861     //----------------------------------------------------------------------
00866     void ioss_add_fields(stk::mesh::Part &part,
00867        stk::mesh::EntityRank part_type,
00868        Ioss::GroupingEntity *entity,
00869        const Ioss::Field::RoleType filter_role,
00870        bool add_all)
00871     {
00872       stk::mesh::MetaData & meta = mesh::MetaData::get(part);
00873       stk::mesh::Part &universal = meta.universal_part();
00874 
00875       const std::vector<mesh::FieldBase*> &fields = meta.get_fields();
00876 
00877       std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
00878       while (I != fields.end()) {
00879   const stk::mesh::FieldBase *f = *I; ++I;
00880   if (stk::io::is_valid_part_field(f, part_type, part, universal, filter_role, add_all)) {
00881     const stk::mesh::FieldBase::Restriction &res = f->restriction(part_type, part);
00882     std::pair<std::string, Ioss::Field::BasicType> field_type;
00883     get_io_field_type(f, res, &field_type);
00884     if (field_type.second != Ioss::Field::INVALID) {
00885       int entity_size = entity->get_property("entity_count").get_int();
00886       const std::string& name = f->name();
00887       entity->field_add(Ioss::Field(name, field_type.second, field_type.first,
00888             filter_role, entity_size));
00889     }
00890   }
00891       }
00892     }
00893 
00903     void define_io_fields(Ioss::GroupingEntity *entity,
00904         Ioss::Field::RoleType role,
00905         stk::mesh::Part &part,
00906         stk::mesh::EntityRank part_type)
00907     {
00908       stk::mesh::MetaData &meta = mesh::MetaData::get(part);
00909 
00910       bool use_cartesian_for_scalar = false;
00911       if (role == Ioss::Field::ATTRIBUTE)
00912   use_cartesian_for_scalar = true;
00913 
00914       Ioss::NameList names;
00915       entity->field_describe(role, &names);
00916 
00917       for (Ioss::NameList::const_iterator I = names.begin(); I != names.end(); ++I) {
00918   // \todo IMPLEMENT Need a field selection mechanism and a field naming
00919   // (ioss_name -> stk::name)  For now, select all and give the
00920   // stk field the same name as the ioss field.
00921 
00922   // Skip the attribute field that is named "attribute"
00923   if (*I == "attribute" && names.size() > 1)
00924     continue;
00925 
00926   // \todo IMPLEMENT Need to determine whether these are
00927   // multi-state fields or constant, or interpolated, or ...
00928   Ioss::Field io_field = entity->get_field(*I);
00929   declare_ioss_field(meta, part_type, part, io_field, use_cartesian_for_scalar);
00930       }
00931     }
00932 
00933     void get_entity_list(Ioss::GroupingEntity *io_entity,
00934        stk::mesh::EntityRank part_type,
00935        const stk::mesh::BulkData &bulk,
00936        std::vector<stk::mesh::Entity*> &entities)
00937     {
00938       std::vector<int> ids ;
00939       io_entity->get_field_data("ids", ids);
00940 
00941       size_t count = ids.size();
00942       entities.resize(count);
00943 
00944       for(size_t i=0; i<count; ++i) {
00945   entities[i] = bulk.get_entity( part_type, ids[i] );
00946       }
00947     }
00948 
00949     void field_data_from_ioss(const stk::mesh::FieldBase *field,
00950             std::vector<stk::mesh::Entity*> &entities,
00951             Ioss::GroupingEntity *io_entity,
00952             const std::string &io_fld_name)
00953     {
00956 
00957       if (field != NULL && io_entity->field_exists(io_fld_name)) {
00958   const Ioss::Field &io_field = io_entity->get_fieldref(io_fld_name);
00959   if (field->type_is<double>()) {
00960     internal_field_data_from_ioss(io_field, field, entities, io_entity,
00961           static_cast<double>(1.0));
00962   } else if (field->type_is<int>()) {
00963     // Make sure the IO field type matches the STK field type.
00964     // By default, all IO fields are created of type 'double'
00965     if (io_field.get_type() != Ioss::Field::INTEGER) {
00966       Ioss::Field &tmp = const_cast<Ioss::Field&>(io_field);
00967       tmp.reset_type(Ioss::Field::INTEGER);
00968     }
00969     internal_field_data_from_ioss(io_field, field, entities, io_entity,
00970           static_cast<int>(1));
00971   }
00972       }
00973     }
00974 
00975     void field_data_to_ioss(const stk::mesh::FieldBase *field,
00976           std::vector<stk::mesh::Entity*> &entities,
00977           Ioss::GroupingEntity *io_entity,
00978           const std::string &io_fld_name,
00979           Ioss::Field::RoleType filter_role)
00980     {
00983 
00984       if (field != NULL && io_entity->field_exists(io_fld_name)) {
00985   const Ioss::Field &io_field = io_entity->get_fieldref(io_fld_name);
00986   if (io_field.get_role() == filter_role) {
00987     if (field->type_is<double>()) {
00988       internal_field_data_to_ioss(io_field, field, entities, io_entity,
00989           static_cast<double>(1.0));
00990     } else if (field->type_is<int>()) {
00991       if (io_field.get_type() != Ioss::Field::INTEGER) {
00992         Ioss::Field &tmp = const_cast<Ioss::Field&>(io_field);
00993         tmp.reset_type(Ioss::Field::INTEGER);
00994       }
00995       internal_field_data_to_ioss(io_field, field, entities, io_entity,
00996           static_cast<int>(1));
00997     }
00998   }
00999       }
01000     }
01001 
01002     //----------------------------------------------------------------------
01003 
01004     // ----------------------------------------------------------------------
01005     // Returns true if 'entity' should be a 'part' in the analysis mesh.
01006     // Returns false if the application is only using a subset of the
01007     // database entities and this entity is not to be used. The
01008     // "omitted" property is set by the application during parsing or
01009     // pre-mesh reading time.
01010     bool include_entity(Ioss::GroupingEntity *entity)
01011     {
01012       assert(entity);
01013 
01014       // Check whether entity has "omitted" property...
01015       bool omitted = (entity->property_exists("omitted")) &&
01016   (entity->get_property("omitted").get_int() == 1);
01017 
01018       return !omitted;
01019     }
01020 
01021     namespace {
01022 
01023       void define_side_block(stk::mesh::Part &part,
01024            Ioss::SideSet *sset,
01025            stk::mesh::EntityRank type,
01026            size_t side_count, int spatial_dimension)
01027       {
01028         const stk::mesh::EntityRank siderank = side_rank(mesh::MetaData::get(part));
01029         const stk::mesh::EntityRank edgerank = edge_rank(mesh::MetaData::get(part));
01030   ThrowRequire(type == siderank || type == edgerank);
01031 
01032   const CellTopologyData *const side_topology =   stk::io::get_cell_topology(part) ?
01033     stk::io::get_cell_topology(part) :
01034     stk::mesh::fem::FEMMetaData::get(part).get_cell_topology(part).getCellTopologyData();
01035 
01036   if (side_topology == NULL) {
01037           std::ostringstream msg ;
01038     msg << " INTERNAL_ERROR: Part " << part.name() << " returned NULL from get_cell_topology()";
01039     throw std::runtime_error( msg.str() );
01040   }
01041 
01042   std::string io_topo = map_topology_cell_to_ioss(side_topology, spatial_dimension);
01043   std::string element_topo_name = "unknown";
01044 
01045   // Get sideblock parent element topology quantities...
01046   // Try to decode from part name...
01047   std::vector<std::string> tokens;
01048         stk::util::tokenize(part.name(), "_", tokens);
01049   if (tokens.size() >= 4) {
01050     // Name of form: "name_eltopo_sidetopo_id" or
01051     //               "name_block_id_sidetopo_id"
01052     // "name" is typically "surface".
01053     const Ioss::ElementTopology *element_topo = Ioss::ElementTopology::factory(tokens[tokens.size()-3], true);
01054     if (element_topo != NULL) {
01055       element_topo_name = element_topo->name();
01056     }
01057   }
01058 
01059   Ioss::SideBlock *side_block = new Ioss::SideBlock( sset->get_database() ,
01060                  part.name() ,
01061                  io_topo, element_topo_name, side_count);
01062   assert(sset->get_side_block(part.name()) == NULL);
01063   sset->add(side_block);
01064 
01065   const mesh::Field<double, mesh::ElementNode> *df = get_distribution_factor_field(part);
01066   if (df != NULL) {
01067     int nodes_per_side = side_topology->node_count;
01068     std::string storage_type = "Real[";
01069     storage_type += Ioss::Utils::to_string(nodes_per_side);
01070     storage_type += "]";
01071     side_block->field_add(Ioss::Field("distribution_factors", Ioss::Field::REAL, storage_type,
01072              Ioss::Field::MESH, side_count));
01073   }
01074       }
01075 
01076       void define_side_blocks(stk::mesh::Part &part,
01077             const stk::mesh::BulkData &bulk_data,
01078             Ioss::SideSet *sset,
01079             stk::mesh::EntityRank type,
01080             int spatial_dimension)
01081       {
01082   mesh::MetaData & meta = mesh::MetaData::get(part);
01083   ThrowRequire(type == face_rank(meta) || type == edge_rank(meta));
01084 
01085   const stk::mesh::PartVector &blocks = part.subsets();
01086   if (blocks.size() > 0) {
01087     for (size_t j = 0; j < blocks.size(); j++) {
01088       mesh::Part & side_block_part = *blocks[j];
01089       stk::mesh::EntityRank side_rank = side_block_part.primary_entity_rank();
01090             mesh::Selector selector = meta.locally_owned_part() & side_block_part;
01091       
01092       size_t num_side = count_selected_entities(selector, bulk_data.buckets(side_rank));
01093 
01094       define_side_block(side_block_part, sset, side_rank, num_side, spatial_dimension);
01095     }
01096   } else {
01097     stk::mesh::EntityRank side_rank = part.primary_entity_rank();
01098           mesh::Selector selector = meta.locally_owned_part() & part;
01099     size_t num_side = count_selected_entities(selector, bulk_data.buckets(side_rank));
01100     define_side_block(part, sset, side_rank, num_side, spatial_dimension);
01101   }
01102       }
01103 
01104       //----------------------------------------------------------------------
01105       //----------------------------------------------------------------------
01106       //----------------------------------------------------------------------
01107 
01108       void define_node_block(stk::mesh::Part &part,
01109            const stk::mesh::BulkData &bulk,
01110            Ioss::Region &io_region)
01111       {
01112   //--------------------------------
01113   // Set the spatial dimension:
01114   mesh::MetaData & meta = mesh::MetaData::get(part);
01115 
01120   mesh::Field<double, mesh::Cartesian> *coord_field =
01121     meta.get_field<stk::mesh::Field<double, mesh::Cartesian> >(std::string("coordinates"));
01122   assert(coord_field != NULL);
01123   const mesh::FieldBase::Restriction &res = coord_field->restriction(node_rank(meta), part);
01124 
01128   const int spatial_dim = res.dimension() ;
01129   io_region.property_add( Ioss::Property("spatial_dimension", spatial_dim));
01130 
01131   //--------------------------------
01132   // Create the special universal node block:
01133 
01134   mesh::Selector selector = meta.locally_owned_part() | meta.globally_shared_part();
01135 
01136   size_t num_nodes = count_selected_entities(selector, bulk.buckets(node_rank(meta)));
01137 
01138   const std::string name("nodeblock_1");
01139 
01140   Ioss::NodeBlock * const nb = new Ioss::NodeBlock(io_region.get_database(),
01141                  name, num_nodes, spatial_dim);
01142   io_region.add( nb );
01143       }
01144 
01145 
01146       void define_element_block(stk::mesh::Part &part,
01147         const stk::mesh::BulkData &bulk,
01148         Ioss::Region &io_region)
01149       {
01150 
01151   mesh::MetaData & meta = mesh::MetaData::get(part);
01152         const stk::mesh::EntityRank elem_rank = element_rank(meta);
01153 
01154         const CellTopologyData * const cell_top =
01155         stk::io::get_cell_topology(part) ?
01156         stk::io::get_cell_topology(part) :
01157         stk::mesh::fem::FEMMetaData::get(part).get_cell_topology(part).getCellTopologyData();
01158 
01159   if (cell_top == NULL) {
01160           std::ostringstream msg ;
01161     msg << " INTERNAL_ERROR: Part " << part.name() << " returned NULL from get_cell_topology()";
01162     throw std::runtime_error( msg.str() );
01163   }
01164 
01165         mesh::Selector selector = meta.locally_owned_part() & part;
01166   const size_t num_elems = count_selected_entities( selector, bulk.buckets(elem_rank));
01167 
01168   int spatial_dim = io_region.get_property("spatial_dimension").get_int();
01169 
01170   // Defer the counting of attributes until after we define the
01171   // element block so we can count them as we add them as fields to
01172   // the element block
01173   Ioss::ElementBlock *eb = new Ioss::ElementBlock(io_region.get_database() ,
01174               part.name() ,
01175               map_topology_cell_to_ioss(cell_top, spatial_dim) ,
01176               num_elems);
01177   io_region.add(eb);
01178 
01179   // Add the attribute fields.
01180   ioss_add_fields(part, part_primary_entity_rank(part), eb, Ioss::Field::ATTRIBUTE);
01181       }
01182 
01183       void define_side_set(stk::mesh::Part &part,
01184          const stk::mesh::BulkData &bulk,
01185          Ioss::Region &io_region)
01186       {
01187         const stk::mesh::EntityRank si_rank = side_rank(mesh::MetaData::get(part));
01188 
01189   bool create_sideset = true;
01190   if (part.subsets().empty()) {
01191     // Only define a sideset for this part if its superset part is
01192     // not a side-containing part..  (i.e., this part is not a subset part
01193     // in a surface...)
01194     const stk::mesh::PartVector &supersets = part.supersets();
01195     for (size_t i=0; i < supersets.size(); i++) {
01196       if (is_part_io_part(*supersets[i]) && supersets[i]->primary_entity_rank() == si_rank) {
01197         create_sideset = false;
01198         break;
01199       }
01200     }
01201   }
01202   if (create_sideset) {
01203     Ioss::SideSet * const ss = new Ioss::SideSet(io_region.get_database(), part.name());
01204 
01205     io_region.add(ss);
01206     int spatial_dim = io_region.get_property("spatial_dimension").get_int();
01207           define_side_blocks(part, bulk, ss, si_rank, spatial_dim);
01208   }
01209       }
01210 
01211       void define_node_set(stk::mesh::Part &part,
01212          const stk::mesh::BulkData &bulk,
01213          Ioss::Region &io_region)
01214       {
01215   mesh::MetaData & meta = mesh::MetaData::get(part);
01216 
01217         mesh::Selector selector = ( meta.locally_owned_part() | meta.globally_shared_part() ) & part;
01218 
01219   const size_t num_nodes =
01220     count_selected_entities(selector, bulk.buckets(node_rank(meta)));
01221 
01222   Ioss::NodeSet * const ns =
01223     new Ioss::NodeSet( io_region.get_database(), part.name(), num_nodes);
01224   io_region.add(ns);
01225       }
01226     } // namespace <blank>
01227 
01228     void define_output_db(Ioss::Region & io_region ,
01229         const mesh::BulkData &bulk_data,
01230         const Ioss::Region *input_region)
01231     {
01232       const mesh::MetaData & meta_data = mesh::MetaData::get(bulk_data);
01233 
01234       const stk::mesh::EntityRank no_rank = node_rank(meta_data);
01235       const stk::mesh::EntityRank el_rank = element_rank(meta_data);
01236       const stk::mesh::EntityRank fa_rank = face_rank(meta_data);
01237       const stk::mesh::EntityRank ed_rank = edge_rank(meta_data);
01238 
01239       io_region.begin_mode( Ioss::STATE_DEFINE_MODEL );
01240 
01241       define_node_block(meta_data.universal_part(), bulk_data, io_region);
01242 
01243       // All parts of the meta data:
01244       const mesh::PartVector & all_parts = meta_data.get_parts();
01245       for (mesh::PartVector::const_iterator i = all_parts.begin();
01246      i != all_parts.end(); ++i) {
01247 
01248   mesh::Part * const part = *i ;
01249 
01250   if (is_part_io_part(*part)) {
01251     if (invalid_rank(part->primary_entity_rank()))
01252       continue;
01253           else if (part->primary_entity_rank() == no_rank)
01254       define_node_set(*part, bulk_data, io_region);
01255           else if (part->primary_entity_rank() == el_rank)
01256       define_element_block(*part, bulk_data, io_region);
01257           else if (part->primary_entity_rank() == fa_rank)
01258       define_side_set(*part, bulk_data, io_region);
01259           else if (part->primary_entity_rank() == ed_rank)
01260       define_side_set(*part, bulk_data, io_region);
01261   }
01262       }
01263 
01264       if (input_region != NULL)
01265   io_region.synchronize_id_and_name(input_region, true);
01266 
01267       io_region.end_mode( Ioss::STATE_DEFINE_MODEL );
01268     }
01269 
01270     //----------------------------------------------------------------------
01271 
01272     namespace {
01273 
01274       size_t get_entities(stk::mesh::Part &part,
01275         const stk::mesh::BulkData &bulk,
01276         std::vector<mesh::Entity*> &entities,
01277         bool include_shared=true)
01278       {
01279   mesh::MetaData & meta = mesh::MetaData::get(part);
01280   mesh::EntityRank type = part_primary_entity_rank(part);
01281   if (invalid_rank(type))
01282     type = node_rank(meta);
01283 
01284   mesh::Selector own_share = meta.locally_owned_part();
01285   if (include_shared)
01286     own_share |= meta.globally_shared_part();
01287   
01288         mesh::Selector selector = part & own_share;
01289   get_selected_entities(selector, bulk.buckets(type), entities);
01290   return entities.size();
01291       }
01292 
01293       void write_side_data_to_ioss( Ioss::GroupingEntity & io ,
01294             mesh::Part * const part ,
01295             const mesh::BulkData & bulk_data )
01296       {
01297         //std::cout << "tmp write_side_data_to_ioss part= " << part->name() << std::endl;
01298   const mesh::MetaData & meta_data = mesh::MetaData::get(*part);
01299 
01300   std::vector<mesh::Entity *> sides ;
01301   size_t num_sides = get_entities(*part, bulk_data, sides, false);
01302 
01303   std::vector<int> side_ids(num_sides);
01304   std::vector<int> elem_side_ids(num_sides*2);
01305 
01306   stk::mesh::EntityRank elem_rank = element_rank(meta_data);
01307   for(size_t i=0; i<num_sides; ++i) {
01308 
01309     const mesh::Entity &side = *sides[i] ;
01310     const mesh::PairIterRelation side_elem = side.relations( elem_rank );
01311 
01312     // Which element to use?
01313     // Any locally owned element that has the "correct" orientation
01314 
01315     const size_t num_side_elem = side_elem.size();
01316 
01317     const mesh::Relation *rel = NULL ;
01318 
01319     for ( size_t j = 0 ; j < num_side_elem && ! rel ; ++j ) {
01320       const mesh::Entity & elem = *side_elem[j].entity();
01321 
01322       if ( elem.bucket().member( meta_data.locally_owned_part() ) &&
01323      (num_side_elem == 1 || stk::mesh::fem::element_side_polarity(elem, side, side_elem[j].identifier())) ) {
01324         rel = &side_elem[j];
01325       }
01326     }
01327 
01328     if (rel == NULL) { // no suitable element found
01329       std::ostringstream oss;
01330       oss << "ERROR, no suitable element found";
01331       throw std::runtime_error(oss.str());
01332     }
01333 
01334     side_ids[i] = side.identifier();
01335 
01336     elem_side_ids[i*2]   = rel->entity()->identifier();
01337     elem_side_ids[i*2+1] = rel->identifier() + 1; // Ioss is 1-based, mesh is 0-based.
01338   }
01339 
01340   const size_t num_ids_written = io.put_field_data("ids", side_ids);
01341   const size_t num_side_written = io.put_field_data("element_side",elem_side_ids);
01342 
01343   if ( num_sides != num_ids_written || num_sides != num_side_written ) {
01344     std::ostringstream msg ;
01345 
01346     msg << "stk::io::write_side_data_to_ioss FAILED for " ;
01347     msg << io.name();
01348     msg << " in Ioss::GroupingEntity::put_field_data:" ;
01349     msg << " num_sides = " << num_sides ;
01350     msg << " , num_ids_written = " << num_ids_written ;
01351     msg << " , num_side_written = " << num_side_written ;
01352     throw std::runtime_error( msg.str() );
01353   }
01354 
01355   const mesh::Field<double, mesh::ElementNode> *df = get_distribution_factor_field(*part);
01356   if (df != NULL) {
01357     field_data_to_ioss(df, sides, &io, "distribution_factors", Ioss::Field::MESH);
01358   }
01359       }
01360 
01361 
01362       //----------------------------------------------------------------------
01363       void output_node_block(Ioss::NodeBlock &nb,
01364            stk::mesh::Part &part,
01365            const stk::mesh::BulkData &bulk)
01366       {
01367   //----------------------------------
01368   // Exactly one node block to obtain the nodal coordinates and ids:
01369   // Note that the "ids" field of the nodes needs to be written
01370   // before any other bulk data that uses node ids since it sets up
01371   // the global->local mapping of nodes for the output database.
01372   // Similarly for the element "ids" field related to bulk data
01373   // using element ids.
01374   std::vector<mesh::Entity *> nodes ;
01375   size_t num_nodes = get_entities(part, bulk, nodes, true);
01376 
01377   std::vector<int> node_ids(num_nodes);
01378   for(size_t i=0; i<num_nodes; ++i) {
01379     const mesh::Entity & node = * nodes[i] ;
01380     node_ids[i] = node.identifier();
01381   }
01382 
01383   size_t num_ids_written = nb.put_field_data("ids", node_ids);
01384   if ( num_nodes != num_ids_written) {
01385     std::ostringstream msg ;
01386     msg << " FAILED in Ioss::NodeBlock::put_field_data:" ;
01387     msg << " num_nodes = " << num_nodes ;
01388     msg << " , num_ids_written = " << num_ids_written ;
01389     throw std::runtime_error( msg.str() );
01390   }
01391 
01396     const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01397     mesh::Field<double, mesh::Cartesian> *coord_field =
01398       meta_data.get_field<stk::mesh::Field<double, mesh::Cartesian> >(std::string("coordinates"));
01399     assert(coord_field != NULL);
01400     field_data_to_ioss(coord_field, nodes, &nb, "mesh_model_coordinates", Ioss::Field::MESH);
01401       }
01402 
01403       void output_element_block(Ioss::ElementBlock *block,
01404         const stk::mesh::BulkData &bulk)
01405       {
01406   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01407   const std::string& name = block->name();
01408   mesh::Part* part = meta_data.get_part(name);
01409 
01410   assert(part != NULL);
01411   std::vector<mesh::Entity *> elements;
01412   size_t num_elems = get_entities(*part, bulk, elements, false);
01413 
01414   const CellTopologyData * cell_topo =
01415               stk::io::get_cell_topology(*part) ?
01416               stk::io::get_cell_topology(*part) :
01417               stk::mesh::fem::FEMMetaData::get(*part).get_cell_topology(*part).getCellTopologyData();
01418   if (cell_topo == NULL) {
01419           std::ostringstream msg ;
01420     msg << " INTERNAL_ERROR: Part " << part->name() << " returned NULL from get_cell_topology()";
01421     throw std::runtime_error( msg.str() );
01422   }
01423   int nodes_per_elem = cell_topo->node_count;
01424 
01425   std::vector<int> elem_ids(num_elems);
01426   std::vector<int> connectivity(num_elems*nodes_per_elem);
01427 
01428   stk::mesh::EntityRank no_rank = node_rank(meta_data);
01429   for (size_t i = 0; i < num_elems; ++i) {
01430 
01431     elem_ids[i] = elements[i]->identifier();
01432 
01433     const mesh::PairIterRelation elem_nodes = elements[i]->relations(no_rank);
01434 
01435     for (int j = 0; j < nodes_per_elem; ++j) {
01436       connectivity[i * nodes_per_elem+j] = elem_nodes[j].entity()->identifier();
01437     }
01438   }
01439 
01440   const size_t num_ids_written = block->put_field_data("ids", elem_ids);
01441   const size_t num_con_written = block->put_field_data("connectivity", connectivity);
01442 
01443   if ( num_elems != num_ids_written || num_elems != num_con_written ) {
01444     std::ostringstream msg ;
01445     msg << " FAILED in Ioss::ElementBlock::put_field_data:" << std::endl ;
01446     msg << "  num_elems = " << num_elems << std::endl ;
01447     msg << "  num_ids_written = " << num_ids_written << std::endl ;
01448     msg << "  num_connectivity_written = " << num_con_written << std::endl ;
01449     throw std::runtime_error( msg.str() );
01450   }
01451 
01452   stk::mesh::EntityRank elem_rank = element_rank(meta_data);
01453   const std::vector<mesh::FieldBase *> &fields = meta_data.get_fields();
01454   std::vector<mesh::FieldBase *>::const_iterator I = fields.begin();
01455   while (I != fields.end()) {
01456     const mesh::FieldBase *f = *I ; ++I ;
01457     const Ioss::Field::RoleType *role = stk::io::get_field_role(*f);
01458     if (role != NULL && *role == Ioss::Field::ATTRIBUTE) {
01459       const mesh::FieldBase::Restriction &res = f->restriction(elem_rank, *part);
01460       if (res.dimension() > 0) {
01461         stk::io::field_data_to_ioss(f, elements, block, f->name(), Ioss::Field::ATTRIBUTE);
01462       }
01463     }
01464   }
01465       }
01466 
01467       void output_node_set(Ioss::NodeSet *ns, const stk::mesh::BulkData &bulk)
01468       {
01469   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01470   const std::string& name = ns->name();
01471   stk::mesh::Part* part = meta_data.get_part(name);
01472   assert(part != NULL);
01473 
01474   std::vector<stk::mesh::Entity *> nodes ;
01475   size_t num_nodes = get_entities(*part, bulk, nodes, true);
01476 
01477   std::vector<int> node_ids(num_nodes);
01478 
01479   for(size_t i=0; i<num_nodes; ++i) {
01480     const stk::mesh::Entity & node = * nodes[i] ;
01481     node_ids[i] = node.identifier();
01482   }
01483   size_t num_ids_written = ns->put_field_data("ids", node_ids);
01484   if ( num_nodes != num_ids_written ) {
01485     std::ostringstream msg ;
01486     msg << " FAILED in Ioss::NodeSet::put_field_data:"
01487         << " num_nodes = " << num_nodes
01488         << ", num_ids_written = " << num_ids_written;
01489     throw std::runtime_error( msg.str() );
01490   }
01491 
01492   stk::mesh::Field<double> *df_field =
01493     meta_data.get_field<stk::mesh::Field<double> >("distribution_factors");
01494   if (df_field != NULL) {
01495     stk::io::field_data_to_ioss(df_field, nodes, ns, "distribution_factors", Ioss::Field::MESH);
01496   }
01497       }
01498 
01499       void output_side_set(Ioss::SideSet *ss,
01500          const stk::mesh::BulkData &bulk)
01501       {
01502   const stk::mesh::MetaData & meta_data = mesh::MetaData::get(bulk);
01503   size_t block_count = ss->block_count();
01504   for (size_t i=0; i < block_count; i++) {
01505     Ioss::SideBlock *block = ss->get_block(i);
01506     if (stk::io::include_entity(block)) {
01507       stk::mesh::Part * const part = meta_data.get_part(block->name());
01508       stk::io::write_side_data_to_ioss(*block, part, bulk);
01509     }
01510   }
01511       }
01512 
01513     } // namespace <blank>
01514 
01515     void write_output_db(Ioss::Region& io_region,
01516        const stk::mesh::BulkData& bulk)
01517     {
01518       const stk::mesh::MetaData & meta = mesh::MetaData::get(bulk);
01519 
01520       io_region.begin_mode( Ioss::STATE_MODEL );
01521 
01522       Ioss::NodeBlock & nb = *io_region.get_node_blocks()[0];
01523       output_node_block(nb, meta.universal_part(), bulk);
01524 
01525       //----------------------------------
01526       const Ioss::ElementBlockContainer& elem_blocks = io_region.get_element_blocks();
01527       for(Ioss::ElementBlockContainer::const_iterator it = elem_blocks.begin();
01528     it != elem_blocks.end(); ++it) {
01529   output_element_block(*it, bulk);
01530       }
01531 
01532       //----------------------------------
01533       const Ioss::NodeSetContainer& node_sets = io_region.get_nodesets();
01534       for(Ioss::NodeSetContainer::const_iterator it = node_sets.begin();
01535     it != node_sets.end(); ++it) {
01536   output_node_set(*it, bulk);
01537       }
01538 
01539       //----------------------------------
01540       const Ioss::SideSetContainer& side_sets = io_region.get_sidesets();
01541       for(Ioss::SideSetContainer::const_iterator it = side_sets.begin();
01542     it != side_sets.end(); ++it) {
01543   output_side_set(*it, bulk);
01544       }
01545 
01546       io_region.end_mode( Ioss::STATE_MODEL );
01547     }
01548 
01549     //----------------------------------------------------------------------
01550     bool is_part_io_part(stk::mesh::Part &part)
01551     {
01552       return NULL != part.attribute<Ioss::GroupingEntity>();
01553     }
01554 
01555     const stk::mesh::Field<double, stk::mesh::ElementNode> *get_distribution_factor_field(const stk::mesh::Part &p)
01556     {
01557       return p.attribute<stk::mesh::Field<double, stk::mesh::ElementNode> >();
01558     }
01559 
01560     void set_distribution_factor_field(stk::mesh::Part &p,
01561                const stk::mesh::Field<double, stk::mesh::ElementNode> &df_field)
01562     {
01563       stk::mesh::MetaData &m = mesh::MetaData::get(p);
01564       m.declare_attribute_no_delete(p,&df_field);
01565     }
01566 
01567     const Ioss::Field::RoleType* get_field_role(const stk::mesh::FieldBase &f)
01568     {
01569       return f.attribute<Ioss::Field::RoleType>();
01570     }
01571 
01572     void set_field_role(stk::mesh::FieldBase &f, const Ioss::Field::RoleType &role)
01573     {
01574       Ioss::Field::RoleType *my_role = new Ioss::Field::RoleType(role);
01575       stk::mesh::MetaData &m = mesh::MetaData::get(f);
01576       const Ioss::Field::RoleType *check = m.declare_attribute_with_delete(f, my_role);
01577       if ( check != my_role ) {
01578         if (*check != *my_role) {
01579           std::ostringstream msg ;
01580     msg << " FAILED in IossBridge -- set_field_role:"
01581         << " The role type had already been set to " << *check
01582         << ", so it is not possible to change it to " << *my_role;
01583     throw std::runtime_error( msg.str() );
01584         }
01585         delete my_role;
01586       }
01587     }
01588 
01589   }//namespace io
01590 }//namespace stk
01591 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends