Sierra Toolkit Version of the Day
FEMHelpers.cpp
00001 #include <stk_mesh/fem/FEMHelpers.hpp>
00002 #include <stk_mesh/fem/FEMMetaData.hpp>
00003 
00004 #include <Shards_CellTopologyTraits.hpp>
00005 
00006 #include <stk_mesh/base/Types.hpp>
00007 #include <stk_mesh/base/BulkData.hpp>
00008 
00009 #include <stk_util/parallel/ParallelReduce.hpp>
00010 
00011 #include <sstream>
00012 #include <stdexcept>
00013 
00014 namespace stk {
00015 namespace mesh {
00016 namespace fem {
00017 
00018 namespace {
00019 
00020 void verify_declare_element_side(
00021     const BulkData & mesh,
00022     const Entity & elem,
00023     const unsigned local_side_id
00024     )
00025 {
00026   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00027 
00028   const CellTopologyData * const side_top =
00029     ( elem_top && local_side_id < elem_top->side_count )
00030     ? elem_top->side[ local_side_id ].topology : NULL ;
00031 
00032   ThrowErrorMsgIf( &mesh != & BulkData::get(elem),
00033     "For elem " << print_entity_key(elem) <<
00034     ", Bulkdata for 'elem' and mesh are different");
00035 
00036   ThrowErrorMsgIf( elem_top && local_side_id >= elem_top->side_count,
00037     "For elem " << print_entity_key(elem) << ", local_side_id " << local_side_id << ", " <<
00038     "local_side_id exceeds " << elem_top->name << ".side_count = " << elem_top->side_count );
00039 
00040   ThrowErrorMsgIf( side_top == NULL,
00041     "For elem " << print_entity_key(elem) << ", local_side_id " << local_side_id << ", " <<
00042     "No element topology found");
00043 }
00044 
00045 } // unnamed namespace
00046 
00047 Entity & declare_element( BulkData & mesh ,
00048                           Part & part ,
00049                           const EntityId elem_id ,
00050                           const EntityId node_id[] )
00051 {
00052   FEMMetaData & fem_meta = FEMMetaData::get(mesh);
00053   const CellTopologyData * const top = fem_meta.get_cell_topology( part ).getCellTopologyData();
00054 
00055   ThrowErrorMsgIf(top == NULL,
00056                   "Part " << part.name() << " does not have a local topology");
00057 
00058   PartVector empty ;
00059   PartVector add( 1 ); add[0] = & part ;
00060 
00061   const EntityRank entity_rank = fem_meta.element_rank();
00062 
00063   Entity & elem = mesh.declare_entity( entity_rank, elem_id, add );
00064 
00065   const EntityRank node_rank = fem_meta.node_rank();
00066 
00067   for ( unsigned i = 0 ; i < top->node_count ; ++i ) {
00068     //declare node if it doesn't already exist
00069     Entity * node = mesh.get_entity( node_rank , node_id[i]);
00070     if ( NULL == node) {
00071       node = & mesh.declare_entity( node_rank , node_id[i], empty );
00072     }
00073 
00074     mesh.declare_relation( elem , *node , i );
00075   }
00076   return elem ;
00077 }
00078 
00079 Entity & declare_element_side(
00080   Entity & elem ,
00081   Entity & side,
00082   const unsigned local_side_id ,
00083   Part * part )
00084 {
00085   BulkData & mesh = BulkData::get(side);
00086 
00087   verify_declare_element_side(mesh, elem, local_side_id);
00088 
00089   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00090 
00091   ThrowErrorMsgIf( elem_top == NULL,
00092       "Element[" << elem.identifier() << "] has no defined topology" );
00093 
00094   const CellTopologyData * const side_top = elem_top->side[ local_side_id ].topology;
00095 
00096   ThrowErrorMsgIf( side_top == NULL,
00097       "Element[" << elem.identifier() << "], local_side_id = " <<
00098       local_side_id << ", side has no defined topology" );
00099 
00100   const unsigned * const side_node_map = elem_top->side[ local_side_id ].node ;
00101 
00102   PartVector add_parts ;
00103 
00104   if ( part ) { add_parts.push_back( part ); }
00105 
00106   mesh.change_entity_parts(side, add_parts);
00107 
00108   mesh.declare_relation( elem , side , local_side_id );
00109 
00110   PairIterRelation rel = elem.relations( FEMMetaData::NODE_RANK );
00111 
00112   for ( unsigned i = 0 ; i < side_top->node_count ; ++i ) {
00113     Entity & node = * rel[ side_node_map[i] ].entity();
00114     mesh.declare_relation( side , node , i );
00115   }
00116 
00117   return side ;
00118 }
00119 
00120 Entity & declare_element_side(
00121   BulkData & mesh ,
00122   const stk::mesh::EntityId global_side_id ,
00123   Entity & elem ,
00124   const unsigned local_side_id ,
00125   Part * part )
00126 {
00127   verify_declare_element_side(mesh, elem, local_side_id);
00128 
00129   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00130 
00131   ThrowErrorMsgIf( elem_top == NULL,
00132       "Element[" << elem.identifier() << "] has no defined topology");
00133 
00134   const CellTopologyData * const side_top = elem_top->side[ local_side_id ].topology;
00135 
00136   ThrowErrorMsgIf( side_top == NULL,
00137       "Element[" << elem.identifier() << "], local_side_id = " <<
00138       local_side_id << ", side has no defined topology" );
00139 
00140   PartVector empty_parts ;
00141 
00142   Entity & side = mesh.declare_entity( side_top->dimension , global_side_id, empty_parts );
00143   return declare_element_side( elem, side, local_side_id, part);
00144 }
00145 
00146 
00147 
00148 const CellTopologyData * get_subcell_nodes(const Entity & entity ,
00149                                            EntityRank subcell_rank ,
00150                                            unsigned subcell_identifier ,
00151                                            EntityVector & subcell_nodes)
00152 {
00153   subcell_nodes.clear();
00154 
00155   // get cell topology
00156   const CellTopologyData* celltopology = get_cell_topology(entity).getCellTopologyData();
00157 
00158   //error checking
00159   {
00160     //no celltopology defined
00161     if (celltopology == NULL) {
00162       return NULL;
00163     }
00164 
00165     // valid ranks fall within the dimension of the cell topology
00166     const bool bad_rank = subcell_rank >= celltopology->dimension;
00167     ThrowInvalidArgMsgIf( bad_rank, "subcell_rank is >= celltopology dimension\n");
00168 
00169     // subcell_identifier must be less than the subcell count
00170     const bool bad_id = subcell_identifier >= celltopology->subcell_count[subcell_rank];
00171     ThrowInvalidArgMsgIf( bad_id,   "subcell_id is >= subcell_count\n");
00172   }
00173 
00174   // Get the cell topology of the subcell
00175   const CellTopologyData * subcell_topology =
00176     celltopology->subcell[subcell_rank][subcell_identifier].topology;
00177 
00178   const int num_nodes_in_subcell = subcell_topology->node_count;
00179 
00180   // For the subcell, get it's local nodes ids
00181   const unsigned* subcell_node_local_ids =
00182     celltopology->subcell[subcell_rank][subcell_identifier].node;
00183 
00184   FEMMetaData & fem_meta = FEMMetaData::get(entity);
00185   const EntityRank node_rank = fem_meta.node_rank();
00186   PairIterRelation node_relations = entity.relations(node_rank);
00187 
00188   subcell_nodes.reserve(num_nodes_in_subcell);
00189 
00190   for (int i = 0; i < num_nodes_in_subcell; ++i ) {
00191     subcell_nodes.push_back( node_relations[subcell_node_local_ids[i]].entity() );
00192   }
00193 
00194   return subcell_topology;
00195 }
00196 
00197 
00198 int get_entity_subcell_id( const Entity & entity ,
00199                            const EntityRank subcell_rank,
00200                            const CellTopologyData * subcell_topology,
00201                            const std::vector<Entity*>& subcell_nodes )
00202 {
00203   const int INVALID_SIDE = -1;
00204 
00205   unsigned num_nodes = subcell_topology->node_count;
00206 
00207   if (num_nodes != subcell_nodes.size()) {
00208     return INVALID_SIDE;
00209   }
00210 
00211   // get topology of elem
00212   const CellTopologyData* entity_topology = get_cell_topology(entity).getCellTopologyData();
00213   if (entity_topology == NULL) {
00214     return INVALID_SIDE;
00215   }
00216 
00217   // get nodal relations for entity
00218   FEMMetaData & fem_meta = FEMMetaData::get(entity);
00219   const EntityRank node_rank = fem_meta.node_rank();
00220   PairIterRelation relations = entity.relations(node_rank);
00221 
00222   const int num_permutations = subcell_topology->permutation_count;
00223 
00224   // Iterate over the subcells of entity...
00225   for (unsigned local_subcell_ordinal = 0;
00226       local_subcell_ordinal < entity_topology->subcell_count[subcell_rank];
00227       ++local_subcell_ordinal) {
00228 
00229     // get topological data for this subcell
00230     const CellTopologyData* curr_subcell_topology =
00231       entity_topology->subcell[subcell_rank][local_subcell_ordinal].topology;
00232 
00233     // If topologies are not the same, there is no way the subcells are the same
00234     if (subcell_topology == curr_subcell_topology) {
00235 
00236       const unsigned* const subcell_node_map = entity_topology->subcell[subcell_rank][local_subcell_ordinal].node;
00237 
00238       // Taking all positive permutations into account, check if this subcell
00239       // has the same nodes as the subcell_nodes argument. Note that this
00240       // implementation preserves the node-order so that we can take
00241       // entity-orientation into account.
00242       for (int p = 0; p < num_permutations; ++p) {
00243 
00244         if (curr_subcell_topology->permutation[p].polarity ==
00245             CELL_PERMUTATION_POLARITY_POSITIVE) {
00246 
00247           const unsigned * const perm_node =
00248             curr_subcell_topology->permutation[p].node ;
00249 
00250           bool all_match = true;
00251           for (unsigned j = 0 ; j < num_nodes; ++j ) {
00252             if (subcell_nodes[j] !=
00253                 relations[subcell_node_map[perm_node[j]]].entity()) {
00254               all_match = false;
00255               break;
00256             }
00257           }
00258 
00259           // all nodes were the same, we have a match
00260           if ( all_match ) {
00261             return local_subcell_ordinal ;
00262           }
00263         }
00264       }
00265     }
00266   }
00267 
00268   return INVALID_SIDE;
00269 }
00270 
00271 bool comm_mesh_counts( BulkData & M ,
00272                        std::vector<size_t> & counts ,
00273                        bool local_flag )
00274 {
00275   const size_t zero = 0 ;
00276 
00277   // Count locally owned entities
00278 
00279   const FEMMetaData & S = FEMMetaData::get(M);
00280   const unsigned entity_rank_count = S.entity_rank_count();
00281   const size_t   comm_count        = entity_rank_count + 1 ;
00282 
00283   std::vector<size_t> local(  comm_count , zero );
00284   std::vector<size_t> global( comm_count , zero );
00285 
00286   ParallelMachine comm = M.parallel();
00287   Part & owns = S.locally_owned_part();
00288 
00289   for ( unsigned i = 0 ; i < entity_rank_count ; ++i ) {
00290     const std::vector<Bucket*> & ks = M.buckets( i );
00291 
00292     std::vector<Bucket*>::const_iterator ik ;
00293 
00294     for ( ik = ks.begin() ; ik != ks.end() ; ++ik ) {
00295       if ( has_superset( **ik , owns ) ) {
00296         local[i] += (*ik)->size();
00297       }
00298     }
00299   }
00300 
00301   local[ entity_rank_count ] = local_flag ;
00302 
00303   stk::all_reduce_sum( comm , & local[0] , & global[0] , comm_count );
00304 
00305   counts.assign( global.begin() , global.begin() + entity_rank_count );
00306 
00307   return 0 < global[ entity_rank_count ] ;
00308 }
00309 
00310 bool element_side_polarity( const Entity & elem ,
00311                             const Entity & side , int local_side_id )
00312 {
00313   // 09/14/10:  TODO:  tscoffe:  Will this work in 1D?
00314   FEMMetaData &fem_meta = FEMMetaData::get(elem);
00315   const bool is_side = side.entity_rank() != fem_meta.edge_rank();
00316   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00317 
00318   const unsigned side_count = ! elem_top ? 0 : (
00319                                 is_side ? elem_top->side_count
00320                                         : elem_top->edge_count );
00321 
00322   ThrowErrorMsgIf( elem_top == NULL,
00323                    "For Element[" << elem.identifier() << "], element has no defined topology");
00324 
00325   ThrowErrorMsgIf( local_side_id < 0 || static_cast<int>(side_count) <= local_side_id,
00326     "For Element[" << elem.identifier() << "], " <<
00327     "side: " << print_entity_key(side) << ", " <<
00328     "local_side_id = " << local_side_id <<
00329     " ; unsupported local_side_id");
00330 
00331   const CellTopologyData * const side_top =
00332     is_side ? elem_top->side[ local_side_id ].topology
00333             : elem_top->edge[ local_side_id ].topology ;
00334 
00335   const unsigned * const side_map =
00336     is_side ? elem_top->side[ local_side_id ].node
00337             : elem_top->edge[ local_side_id ].node ;
00338 
00339   const PairIterRelation elem_nodes = elem.relations( FEMMetaData::NODE_RANK );
00340   const PairIterRelation side_nodes = side.relations( FEMMetaData::NODE_RANK );
00341 
00342   bool good = true ;
00343   for ( unsigned j = 0 ; good && j < side_top->node_count ; ++j ) {
00344     good = side_nodes[j].entity() == elem_nodes[ side_map[j] ].entity();
00345   }
00346   return good ;
00347 }
00348 
00349 }
00350 }
00351 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends