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<11> >()), m_element_rank);
00154     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Tetrahedron<8> >()), m_element_rank);
00155 
00156     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<5> >()), m_element_rank);
00157     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<13> >()), m_element_rank);
00158     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Pyramid<14> >()), m_element_rank);
00159 
00160     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<6> >()), m_element_rank);
00161     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<15> >()), m_element_rank);
00162     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Wedge<18> >()), m_element_rank);
00163 
00164     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<8> >()), m_element_rank);
00165     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<20> >()), m_element_rank);
00166     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::Hexahedron<27> >()), m_element_rank);
00167 
00168     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellTriangle<3> >()), m_element_rank);
00169     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellTriangle<6> >()), m_element_rank);
00170 
00171     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<4> >()), m_element_rank);
00172     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<8> >()), m_element_rank);
00173     register_cell_topology(fem::CellTopology(shards::getCellTopologyData< shards::ShellQuadrilateral<9> >()), m_element_rank);
00174   }
00175 }
00176 
00177 void FEMMetaData::register_cell_topology(const fem::CellTopology cell_topology, EntityRank entity_rank)
00178 {
00179   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::register_cell_topology: FEM_initialize() must be called before this function");
00180 
00181   CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
00182 
00183   const bool       duplicate     = it != m_cellTopologyPartEntityRankMap.end();
00184   const EntityRank existing_rank = duplicate ? (*it).second.second : 0;
00185 
00186   ThrowInvalidArgMsgIf(m_spatial_dimension < entity_rank,
00187     "entity_rank " << entity_rank << ", " <<
00188     "exceeds maximum spatial_dimension = " << m_spatial_dimension );
00189 
00190   ThrowErrorMsgIf(duplicate && existing_rank != entity_rank,
00191     "For args: cell_topolgy " << cell_topology.getName() << " and entity_rank " << entity_rank << ", " <<
00192     "previously declared rank = " << existing_rank );
00193 
00194   if (! duplicate) {
00195     std::string part_name = std::string("{FEM_ROOT_CELL_TOPOLOGY_PART_") + std::string(cell_topology.getName()) + std::string("}");
00196 
00197     ThrowErrorMsgIf(get_part(part_name) != 0, "Cannot register topology with same name as existing part '" << cell_topology.getName() << "'" );
00198 
00199     Part &part = declare_part(part_name, entity_rank);
00200     m_cellTopologyPartEntityRankMap[cell_topology] = CellTopologyPartEntityRankMap::mapped_type(&part, entity_rank);
00201 
00202     assign_cell_topology(m_partCellTopologyVector, part.mesh_meta_data_ordinal(), cell_topology);
00203   }
00204   //check_topo_db();
00205 }
00206 
00207 
00208 fem::CellTopology
00209 FEMMetaData::get_cell_topology(
00210   const std::string &   topology_name) const
00211 {
00212   std::string part_name = std::string("{FEM_ROOT_CELL_TOPOLOGY_PART_") + topology_name + std::string("}");
00213 
00214   Part *part = get_part(part_name);
00215   if (part)
00216     return get_cell_topology(*part);
00217   else
00218     return fem::CellTopology();
00219 }
00220 
00221 
00222 Part &FEMMetaData::get_cell_topology_root_part(const fem::CellTopology cell_topology) const
00223 {
00224   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::get_cell_topology_root_part: FEM_initialize() must be called before this function");
00225   CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
00226   ThrowErrorMsgIf(it == m_cellTopologyPartEntityRankMap.end(),
00227                   "Cell topology " << cell_topology.getName() <<
00228                   " has not been registered");
00229 
00230   return *(*it).second.first;
00231 }
00232 
00238 fem::CellTopology FEMMetaData::get_cell_topology( const Part & part) const
00239 {
00240   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::get_cell_topology: FEM_initialize() must be called before this function");
00241   fem::CellTopology cell_topology;
00242 
00243   PartOrdinal part_ordinal = part.mesh_meta_data_ordinal();
00244   if (part_ordinal < m_partCellTopologyVector.size())
00245     {
00246       cell_topology = m_partCellTopologyVector[part_ordinal];
00247     }
00248 
00249   return cell_topology;
00250 }
00251 
00252 #if 0
00253   void FEMMetaData::check_topo_db()
00254   {
00255     std::cout << "FEMMetaData::check_topo_db... m_partCellTopologyVector.size() = " << m_partCellTopologyVector.size() <<  std::endl;
00256 
00257   fem::CellTopology cell_topology;
00258 
00259   for (unsigned i = 0; i <  m_partCellTopologyVector.size(); i++)
00260     {
00261       cell_topology = m_partCellTopologyVector[i];
00262   if (!cell_topology.getCellTopologyData())
00263     {
00264       std::cout << "bad topology in FEMMetaData::check_topo_db" << std::endl;
00265     }
00266       ThrowRequireMsg(cell_topology.getCellTopologyData(), "bad topology in FEMMetaData::check_topo_db");
00267 
00268     }
00269     std::cout << "FEMMetaData::check_topo_db...done" << std::endl;
00270 
00271   }
00272 #endif
00273 
00274 namespace {
00275 
00276 bool root_part_in_subset(stk::mesh::Part & part)
00277 {
00278   if (is_cell_topology_root_part(part)) {
00279     return true;
00280   }
00281   const PartVector & subsets = part.subsets();
00282   for (PartVector::const_iterator it=subsets.begin() ; it != subsets.end() ; ++it) {
00283     if (is_cell_topology_root_part( **it )) {
00284       return true;
00285     }
00286   }
00287   return false;
00288 }
00289 
00290 void find_cell_topologies_in_part_and_subsets_of_same_rank(const Part & part, EntityRank rank, std::set<fem::CellTopology> & topologies_found)
00291 {
00292   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(part);
00293   fem::CellTopology top = fem_meta.get_cell_topology(part);
00294   if ((top.isValid() && (part.primary_entity_rank() == rank))) {
00295     topologies_found.insert(top);
00296   }
00297   const PartVector & subsets = part.subsets();
00298   for (PartVector::const_iterator it=subsets.begin() ; it != subsets.end() ; ++it) {
00299     top = fem_meta.get_cell_topology(**it);
00300     if (top.isValid() && ( (**it).primary_entity_rank() == rank) ) {
00301       topologies_found.insert(top);
00302     }
00303   }
00304 }
00305 
00306 } // namespace
00307 
00308 void FEMMetaData::declare_part_subset( Part & superset , Part & subset )
00309 {
00310   ThrowRequireMsg(is_FEM_initialized(),"FEMMetaData::declare_part_subset: FEM_initialize() must be called before this function");
00311   fem::CellTopology superset_top = get_cell_topology(superset);
00312 
00313   const bool no_superset_topology = !superset_top.isValid();
00314   if ( no_superset_topology ) {
00315     m_meta_data.declare_part_subset(superset,subset);
00316     return;
00317   }
00318   // Check for cell topology root parts in subset or subset's subsets
00319   const bool subset_has_root_part = root_part_in_subset(subset);
00320   ThrowErrorMsgIf( subset_has_root_part, "FEMMetaData::declare_part_subset:  Error, root cell topology part found in subset or below." );
00321 
00322   std::set<fem::CellTopology> cell_topologies;
00323   find_cell_topologies_in_part_and_subsets_of_same_rank(subset,superset.primary_entity_rank(),cell_topologies);
00324 
00325   ThrowErrorMsgIf( cell_topologies.size() > 1,
00326       "FEMMetaData::declare_part_subset:  Error, multiple cell topologies of rank "
00327       << superset.primary_entity_rank()
00328       << " defined below subset"
00329       );
00330   const bool non_matching_cell_topology = ((cell_topologies.size() == 1) && (*cell_topologies.begin() != superset_top));
00331   ThrowErrorMsgIf( non_matching_cell_topology,
00332       "FEMMetaData::declare_part_subset:  Error, superset topology = "
00333       << superset_top.getName() << " does not match the topology = "
00334       << cell_topologies.begin()->getName()
00335       << " coming from the subset part"
00336       );
00337   // Everything is Okay!
00338   m_meta_data.declare_part_subset(superset,subset);
00339   // Update PartCellTopologyVector for "subset" and same-rank subsets, ad nauseum
00340   if (subset.primary_entity_rank() == superset.primary_entity_rank()) {
00341     assign_cell_topology(m_partCellTopologyVector, subset.mesh_meta_data_ordinal(), superset_top);
00342     const PartVector & subset_parts = subset.subsets();
00343     for (PartVector::const_iterator it=subset_parts.begin() ; it != subset_parts.end() ; ++it) {
00344       const Part & it_part = **it;
00345       if (it_part.primary_entity_rank() == superset.primary_entity_rank()) {
00346         assign_cell_topology(m_partCellTopologyVector, it_part.mesh_meta_data_ordinal(), superset_top);
00347       }
00348     }
00349   }
00350 }
00351 
00352 EntityRank FEMMetaData::get_entity_rank(
00353   const fem::CellTopology       cell_topology) const
00354 {
00355   CellTopologyPartEntityRankMap::const_iterator it = m_cellTopologyPartEntityRankMap.find(cell_topology);
00356   if (it == m_cellTopologyPartEntityRankMap.end())
00357     return INVALID_RANK;
00358   else
00359     return (*it).second.second;
00360 }
00361 
00362 //--------------------------------------------------------------------------------
00363 // Free Functions
00364 //--------------------------------------------------------------------------------
00365 
00366 bool is_cell_topology_root_part(const Part & part) {
00367   fem::FEMMetaData & fem_meta = fem::FEMMetaData::get(part);
00368   fem::CellTopology top = fem_meta.get_cell_topology(part);
00369   if (top.isValid()) {
00370     const Part & root_part = fem_meta.get_cell_topology_root_part(top);
00371     return (root_part == part);
00372   }
00373   return false;
00374 }
00375 
00381 void set_cell_topology(
00382   Part &                        part,
00383   fem::CellTopology             cell_topology)
00384 {
00385   FEMMetaData& fem_meta = FEMMetaData::get(part);
00386 
00387   ThrowRequireMsg(fem_meta.is_FEM_initialized(),"set_cell_topology: FEM_initialize() must be called before this function");
00388 
00389   Part &root_part = fem_meta.get_cell_topology_root_part(cell_topology);
00390   fem_meta.declare_part_subset(root_part, part);
00391 }
00392 
00393 std::vector<std::string>
00394 entity_rank_names( size_t spatial_dimension )
00395 {
00396   ThrowInvalidArgMsgIf( spatial_dimension < 1 || 3 < spatial_dimension,
00397                         "Invalid spatial dimension = " << spatial_dimension );
00398 
00399   std::vector< std::string > names ;
00400 
00401   names.reserve( spatial_dimension + 1 );
00402 
00403   names.push_back( std::string( "NODE" ) );
00404 
00405   if ( 1 < spatial_dimension ) { names.push_back( std::string("EDGE") ); }
00406   if ( 2 < spatial_dimension ) { names.push_back( std::string("FACE") ); }
00407 
00408   names.push_back( std::string("ELEMENT") );
00409 
00410   return names ;
00411 }
00412 
00413 
00414 CellTopology
00415 get_cell_topology(
00416   const Bucket &                bucket)
00417 {
00418   const BulkData   & bulk_data = BulkData::get(bucket);
00419   const MetaData   & meta_data = MetaData::get(bulk_data);
00420   const PartVector & all_parts = meta_data.get_parts();
00421 
00422   FEMMetaData &fem = FEMMetaData::get(meta_data);
00423 
00424   CellTopology cell_topology;
00425 
00426   const std::pair< const unsigned *, const unsigned * > supersets = bucket.superset_part_ordinals();
00427 
00428   if (supersets.first != supersets.second) {
00429     const Part *first_found_part = 0;
00430 
00431     for ( const unsigned * it = supersets.first ; it != supersets.second ; ++it ) {
00432 
00433       const Part & part = * all_parts[*it] ;
00434 
00435       if ( part.primary_entity_rank() == bucket.entity_rank() ) {
00436 
00437         CellTopology top = fem.get_cell_topology( part );
00438 
00439         if ( ! cell_topology.getCellTopologyData() ) {
00440           cell_topology = top ;
00441 
00442           if (!first_found_part)
00443             first_found_part = &part;
00444         }
00445         else {
00446           ThrowErrorMsgIf( top.getCellTopologyData() && top != cell_topology,
00447             "Cell topology is ambiguously defined. It is defined as " << cell_topology.getName() <<
00448             " on part " << first_found_part->name() << " and as " << top.getName() << " on its superset part " << part.name() );
00449         }
00450       }
00451     }
00452   }
00453 
00454   return cell_topology ;
00455 }
00456 
00457 } // namespace fem
00458 } // namespace mesh
00459 } // namespace stk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends