Sierra Toolkit Version of the Day
FEMMetaData.cpp
00001 #include <set>
00002 #include <stk_mesh/fem/FEMMetaData.hpp>
00003 
00004 #include <Shards_CellTopology.hpp>
00005 #include <stk_mesh/base/BulkData.hpp>
00006 #include <stk_mesh/base/Bucket.hpp>
00007 #include <stk_mesh/base/Entity.hpp>
00008 #include <stk_mesh/base/Ghosting.hpp>
00009 
00010 namespace stk {
00011 namespace mesh {
00012 namespace fem {
00013 
00014 namespace {
00015 
00016 void assign_cell_topology(
00017   FEMMetaData::PartCellTopologyVector & part_cell_topology_vector,
00018   size_t                   part_ordinal,
00019   const fem::CellTopology  cell_topology)
00020 {
00021   if (part_ordinal >= part_cell_topology_vector.size())
00022     part_cell_topology_vector.resize(part_ordinal + 1);
00023 
00024   part_cell_topology_vector[part_ordinal] = cell_topology;
00025 
00026   if (!cell_topology.getCellTopologyData())
00027     {
00028       std::cout << "bad topology in FEMMetaData::assign_cell_topology" << std::endl;
00029     }
00030 
00031   ThrowRequireMsg(cell_topology.getCellTopologyData(), "bad topology in FEMMetaData::assign_cell_topology");
00032 }
00033 
00034 } // namespace
00035 
00036 FEMMetaData::FEMMetaData()
00037 :
00038     m_fem_initialized(false),
00039     m_spatial_dimension(0),
00040     m_side_rank(INVALID_RANK),
00041     m_element_rank(INVALID_RANK)
00042 {
00043   // Attach FEMMetaData as attribute on MetaData to enable "get accessors" to FEMMetaData
00044   m_meta_data.declare_attribute_no_delete<FEMMetaData>(this);
00045 }
00046 
00047 FEMMetaData::FEMMetaData(size_t spatial_dimension,
00048                          const std::vector<std::string>& in_entity_rank_names)
00049   :
00050     m_fem_initialized(false),
00051     m_spatial_dimension(0),
00052     m_side_rank(INVALID_RANK),
00053     m_element_rank(INVALID_RANK)
00054 {
00055   // Attach FEMMetaData as attribute on MetaData to enable "get accessors" to FEMMetaData
00056   m_meta_data.declare_attribute_no_delete<FEMMetaData>(this);
00057 
00058   FEM_initialize(spatial_dimension, in_entity_rank_names);
00059 }
00060 
00061 void FEMMetaData::FEM_initialize(size_t spatial_dimension, const std::vector<std::string>& rank_names)
00062 {
00063   ThrowRequireMsg(!m_fem_initialized,"FEM functionality in FEMMetaData can only be initialized once.");
00064   if ( rank_names.empty() ) {
00065     m_entity_rank_names = fem::entity_rank_names(spatial_dimension);
00066   }
00067   else {
00068     ThrowRequireMsg(rank_names.size() >= spatial_dimension+1,
00069                     "Entity rank name vector must name every rank");
00070     m_entity_rank_names = rank_names;
00071   }
00072   internal_set_spatial_dimension_and_ranks(spatial_dimension);
00073   m_meta_data.set_entity_rank_names(m_entity_rank_names);
00074   m_fem_initialized = true;
00075   internal_declare_known_cell_topology_parts();
00076 }
00077 
00078 void FEMMetaData::internal_set_spatial_dimension_and_ranks(size_t spatial_dimension)
00079 {
00080   ThrowRequireMsg( spatial_dimension != 0, "FEMMetaData::internal_set_spatial_dimension_and_ranks: spatial_dimension == 0!");
00081   m_spatial_dimension = spatial_dimension;
00082 
00083   // TODO:  Decide on correct terminology for the FEM Entity Ranks (consider topological vs spatial names and incompatibilities).
00084   // spatial_dimension = 1
00085   // node = 0, edge = 0, face = 0, side = 0, element = 1
00086   // spatial_dimension = 2
00087   // node = 0, edge = 1, face = 1, side = 1, element = 2
00088   // spatial_dimension = 3
00089   // node = 0, edge = 1, face = 2, side = 2, element = 3
00090   // spatial_dimension = 4
00091   // node = 0, edge = 1, face = 2, side = 3, element = 4
00092   m_side_rank = m_spatial_dimension - 1;
00093   m_element_rank = m_spatial_dimension;
00094 
00095 }
00096 
00097 void FEMMetaData::internal_declare_known_cell_topology_parts()
00098 {
00099   // Load up appropriate standard cell topologies.
00100   register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Node >()), NODE_RANK);
00101 
00102   if (m_spatial_dimension == 1) {
00103 
00104     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Particle >()), m_element_rank);
00105 
00106     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<2> >()), m_element_rank); // ???
00107     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<3> >()), m_element_rank); // ???
00108 
00109   }
00110 
00111   else if (m_spatial_dimension == 2) {
00112 
00113     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<2> >()), m_side_rank);
00114     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<3> >()), m_side_rank);
00115 
00116     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Particle >()), m_element_rank);
00117 
00118     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<3> >()), m_element_rank);
00119     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<6> >()), m_element_rank);
00120     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<4> >()), m_element_rank);
00121 
00122     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()), m_element_rank);
00123     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<8> >()), m_element_rank);
00124     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<9> >()), m_element_rank);
00125 
00126     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<2> >()), m_element_rank);
00127     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<3> >()), m_element_rank);
00128 
00129     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellLine<2> >()), m_element_rank);
00130     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellLine<3> >()), m_element_rank);
00131   }
00132 
00133   else if (m_spatial_dimension == 3) {
00134 
00135     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<2> >()), EDGE_RANK);
00136     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Line<3> >()), EDGE_RANK);
00137 
00138     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<3> >()), m_side_rank);
00139     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<6> >()), m_side_rank);
00140     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Triangle<4> >()), m_side_rank);
00141 
00142     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()), m_side_rank);
00143     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<8> >()), m_side_rank);
00144     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<9> >()), m_side_rank);
00145 
00146     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Particle >()), m_element_rank);
00147 
00148     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<2> >()), m_element_rank);
00149     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Beam<3> >()), m_element_rank);
00150 
00151     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<4> >()), m_element_rank);
00152     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<10> >()), m_element_rank);
00153     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<8> >()), m_element_rank);
00154 
00155     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<5> >()), m_element_rank);
00156     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<13> >()), m_element_rank);
00157     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<14> >()), m_element_rank);
00158 
00159     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<6> >()), m_element_rank);
00160     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<15> >()), m_element_rank);
00161     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<18> >()), m_element_rank);
00162 
00163     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<8> >()), m_element_rank);
00164     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<20> >()), m_element_rank);
00165     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<27> >()), m_element_rank);
00166 
00167     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellTriangle<3> >()), m_element_rank);
00168     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellTriangle<6> >()), m_element_rank);
00169 
00170     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<4> >()), m_element_rank);
00171     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<8> >()), m_element_rank);
00172     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<9> >()), m_element_rank);
00173   }
00174 }
00175 
00176 void FEMMetaData::register_cell_topology(const fem::CellTopology cell_topology, EntityRank entity_rank)
00177 {
00178   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::register_cell_topology: FEM_initialize() must be called before this function");
00179 
00180   CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
00181 
00182   const bool       duplicate     = it != m_cellTopologyPartEntityRankMap.end();
00183   const EntityRank existing_rank = duplicate ? (*it).second.second : 0;
00184 
00185   ThrowInvalidArgMsgIf(m_spatial_dimension < entity_rank,
00186     "entity_rank " << entity_rank << ", " <<
00187     "exceeds maximum spatial_dimension = " << m_spatial_dimension );
00188 
00189   ThrowErrorMsgIf(duplicate && existing_rank != entity_rank,
00190     "For args: cell_topolgy " << cell_topology.getName() << " and entity_rank " << entity_rank << ", " <<
00191     "previously declared rank = " << existing_rank );
00192 
00193   if (! duplicate) {
00194     std::string part_name = std::string("{FEM_ROOT_CELL_TOPOLOGY_PART_") + std::string(cell_topology.getName()) + std::string("}");
00195 
00196     ThrowErrorMsgIf(get_part(part_name) != 0, "Cannot register topology with same name as existing part '" << cell_topology.getName() << "'" );
00197 
00198     Part &part = declare_part(part_name, entity_rank);
00199     m_cellTopologyPartEntityRankMap[cell_topology] = CellTopologyPartEntityRankMap::mapped_type(&part, entity_rank);
00200 
00201     assign_cell_topology(m_partCellTopologyVector, part.mesh_meta_data_ordinal(), cell_topology);
00202   }
00203   //check_topo_db();
00204 }
00205 
00206 
00207 fem::CellTopology
00208 FEMMetaData::get_cell_topology(
00209   const std::string &   topology_name) const
00210 {
00211   std::string part_name = std::string("{FEM_ROOT_CELL_TOPOLOGY_PART_") + topology_name + std::string("}");
00212 
00213   Part *part = get_part(part_name);
00214   if (part)
00215     return get_cell_topology(*part);
00216   else
00217     return fem::CellTopology();
00218 }
00219 
00220 
00221 Part &FEMMetaData::get_cell_topology_root_part(const fem::CellTopology cell_topology) const
00222 {
00223   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::get_cell_topology_root_part: FEM_initialize() must be called before this function");
00224   CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
00225   ThrowErrorMsgIf(it == m_cellTopologyPartEntityRankMap.end(),
00226                   "Cell topology " << cell_topology.getName() <<
00227                   " has not been registered");
00228 
00229   return *(*it).second.first;
00230 }
00231 
00237 fem::CellTopology FEMMetaData::get_cell_topology( const Part & part) const
00238 {
00239   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::get_cell_topology: FEM_initialize() must be called before this function");
00240   fem::CellTopology cell_topology;
00241 
00242   PartOrdinal part_ordinal = part.mesh_meta_data_ordinal();
00243   if (part_ordinal < m_partCellTopologyVector.size())
00244     {
00245       cell_topology = m_partCellTopologyVector[part_ordinal];
00246     }
00247 
00248   return cell_topology;
00249 }
00250 
00251 #if 0
00252   void FEMMetaData::check_topo_db()
00253   {
00254     std::cout << "FEMMetaData::check_topo_db... m_partCellTopologyVector.size() = " << m_partCellTopologyVector.size() <<  std::endl;
00255 
00256   fem::CellTopology cell_topology;
00257 
00258   for (unsigned i = 0; i <  m_partCellTopologyVector.size(); i++)
00259     {
00260       cell_topology = m_partCellTopologyVector[i];
00261   if (!cell_topology.getCellTopologyData())
00262     {
00263       std::cout << "bad topology in FEMMetaData::check_topo_db" << std::endl;
00264     }
00265       ThrowRequireMsg(cell_topology.getCellTopologyData(), "bad topology in FEMMetaData::check_topo_db");
00266 
00267     }
00268     std::cout << "FEMMetaData::check_topo_db...done" << std::endl;
00269 
00270   }
00271 #endif
00272 
00273 namespace {
00274 
00275 bool root_part_in_subset(stk::mesh::Part & part)
00276 {
00277   if (is_cell_topology_root_part(part)) {
00278     return true;
00279   }
00280   const PartVector & subsets = part.subsets();
00281   for (PartVector::const_iterator it=subsets.begin() ; it != subsets.end() ; ++it) {
00282     if (is_cell_topology_root_part( **it )) {
00283       return true;
00284     }
00285   }
00286   return false;
00287 }
00288 
00289 void find_cell_topologies_in_part_and_subsets_of_same_rank(const Part & part, EntityRank rank, std::set<fem::CellTopology> & topologies_found)
00290 {
00291   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(part);
00292   fem::CellTopology top = fem_meta.get_cell_topology(part);
00293   if ((top.isValid() && (part.primary_entity_rank() == rank))) {
00294     topologies_found.insert(top);
00295   }
00296   const PartVector & subsets = part.subsets();
00297   for (PartVector::const_iterator it=subsets.begin() ; it != subsets.end() ; ++it) {
00298     top = fem_meta.get_cell_topology(**it);
00299     if (top.isValid() && ( (**it).primary_entity_rank() == rank) ) {
00300       topologies_found.insert(top);
00301     }
00302   }
00303 }
00304 
00305 } // namespace
00306 
00307 void FEMMetaData::declare_part_subset( Part & superset , Part & subset )
00308 {
00309   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::declare_part_subset: FEM_initialize() must be called before this function");
00310   fem::CellTopology superset_top = get_cell_topology(superset);
00311 
00312   const bool no_superset_topology = !superset_top.isValid();
00313   if ( no_superset_topology ) {
00314     m_meta_data.declare_part_subset(superset,subset);
00315     return;
00316   }
00317   // Check for cell topology root parts in subset or subset's subsets
00318   const bool subset_has_root_part = root_part_in_subset(subset);
00319   ThrowErrorMsgIf( subset_has_root_part, "FEMMetaData::declare_part_subset:  Error, root cell topology part found in subset or below." );
00320 
00321   std::set<fem::CellTopology> cell_topologies;
00322   find_cell_topologies_in_part_and_subsets_of_same_rank(subset,superset.primary_entity_rank(),cell_topologies);
00323 
00324   ThrowErrorMsgIf( cell_topologies.size() > 1,
00325       "FEMMetaData::declare_part_subset:  Error, multiple cell topologies of rank "
00326       << superset.primary_entity_rank()
00327       << " defined below subset"
00328       );
00329   const bool non_matching_cell_topology = ((cell_topologies.size() == 1) && (*cell_topologies.begin() != superset_top));
00330   ThrowErrorMsgIf( non_matching_cell_topology,
00331       "FEMMetaData::declare_part_subset:  Error, superset topology = "
00332       << superset_top.getName() << " does not match the topology = "
00333       << cell_topologies.begin()->getName()
00334       << " coming from the subset part"
00335       );
00336   // Everything is Okay!
00337   m_meta_data.declare_part_subset(superset,subset);
00338   // Update PartCellTopologyVector for "subset" and same-rank subsets, ad nauseum
00339   if (subset.primary_entity_rank() == superset.primary_entity_rank()) {
00340     assign_cell_topology(m_partCellTopologyVector, subset.mesh_meta_data_ordinal(), superset_top);
00341     const PartVector & subset_parts = subset.subsets();
00342     for (PartVector::const_iterator it=subset_parts.begin() ; it != subset_parts.end() ; ++it) {
00343       const Part & it_part = **it;
00344       if (it_part.primary_entity_rank() == superset.primary_entity_rank()) {
00345         assign_cell_topology(m_partCellTopologyVector, it_part.mesh_meta_data_ordinal(), superset_top);
00346       }
00347     }
00348   }
00349 }
00350 
00351 EntityRank FEMMetaData::get_entity_rank(
00352   const fem::CellTopology       cell_topology) const
00353 {
00354   CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
00355   if (it == m_cellTopologyPartEntityRankMap.end())
00356     return INVALID_RANK;
00357   else
00358     return (*it).second.second;
00359 }
00360 
00361 //--------------------------------------------------------------------------------
00362 // Free Functions
00363 //--------------------------------------------------------------------------------
00364 
00365 bool is_cell_topology_root_part(const Part & part) {
00366   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(part);
00367   fem::CellTopology top = fem_meta.get_cell_topology(part);
00368   if (top.isValid()) {
00369     const Part & root_part = fem_meta.get_cell_topology_root_part(top);
00370     return (root_part == part);
00371   }
00372   return false;
00373 }
00374 
00380 void set_cell_topology(
00381   Part &                        part,
00382   fem::CellTopology             cell_topology)
00383 {
00384   FEMMetaData& fem_meta = FEMMetaData::get(part);
00385 
00386   ThrowRequireMsg(fem_meta.is_FEM_initialized(),"set_cell_topology: FEM_initialize() must be called before this function");
00387 
00388   Part &root_part = fem_meta.get_cell_topology_root_part(cell_topology);
00389   fem_meta.declare_part_subset(root_part, part);
00390 }
00391 
00392 std::vector<std::string>
00393 entity_rank_names( size_t spatial_dimension )
00394 {
00395   ThrowInvalidArgMsgIf( spatial_dimension < 1 || 3 < spatial_dimension,
00396                         "Invalid spatial dimension = " << spatial_dimension );
00397 
00398   std::vector< std::string > names ;
00399 
00400   names.reserve( spatial_dimension + 1 );
00401 
00402   names.push_back( std::string( "NODE" ) );
00403 
00404   if ( 1 < spatial_dimension ) { names.push_back( std::string("EDGE") ); }
00405   if ( 2 < spatial_dimension ) { names.push_back( std::string("FACE") ); }
00406 
00407   names.push_back( std::string("ELEMENT") );
00408 
00409   return names ;
00410 }
00411 
00412 
00413 CellTopology
00414 get_cell_topology(
00415   const Bucket &                bucket)
00416 {
00417   const BulkData   & bulk_data = BulkData::get(bucket);
00418   const MetaData   & meta_data = MetaData::get(bulk_data);
00419   const PartVector & all_parts = meta_data.get_parts();
00420 
00421   FEMMetaData &fem = FEMMetaData::get(meta_data);
00422 
00423   CellTopology cell_topology;
00424 
00425   const std::pair< const unsigned *, const unsigned * > supersets = bucket.superset_part_ordinals();
00426 
00427   if (supersets.first != supersets.second) {
00428     const Part *first_found_part = 0;
00429 
00430     for ( const unsigned * it = supersets.first ; it != supersets.second ; ++it ) {
00431 
00432       const Part & part = * all_parts[*it] ;
00433 
00434       if ( part.primary_entity_rank() == bucket.entity_rank() ) {
00435 
00436         CellTopology top = fem.get_cell_topology( part );
00437 
00438         if ( ! cell_topology.getCellTopologyData() ) {
00439           cell_topology = top ;
00440 
00441           if (!first_found_part)
00442             first_found_part = &part;
00443         }
00444         else {
00445           ThrowErrorMsgIf( top.getCellTopologyData() && top != cell_topology,
00446             "Cell topology is ambiguously defined. It is defined as " << cell_topology.getName() <<
00447             " on part " << first_found_part->name() << " and as " << top.getName() << " on its superset part " << part.name() );
00448         }
00449       }
00450     }
00451   }
00452 
00453   return cell_topology ;
00454 }
00455 
00456 } // namespace fem
00457 } // namespace mesh
00458 } // namespace stk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends