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