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 void verify_declare_element_edge(
00046     const BulkData & mesh,
00047     const Entity & elem,
00048     const unsigned local_edge_id
00049     )
00050 {
00051   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00052 
00053   const CellTopologyData * const edge_top =
00054     ( elem_top && local_edge_id < elem_top->edge_count )
00055     ? elem_top->edge[ local_edge_id ].topology : NULL ;
00056 
00057   ThrowErrorMsgIf( &mesh != & BulkData::get(elem),
00058     "For elem " << print_entity_key(elem) <<
00059     ", Bulkdata for 'elem' and mesh are different");
00060 
00061   ThrowErrorMsgIf( elem_top && local_edge_id >= elem_top->edge_count,
00062     "For elem " << print_entity_key(elem) << ", local_edge_id " << local_edge_id << ", " <<
00063     "local_edge_id exceeds " << elem_top->name << ".edge_count = " << elem_top->edge_count );
00064 
00065   ThrowErrorMsgIf( edge_top == NULL,
00066     "For elem " << print_entity_key(elem) << ", local_edge_id " << local_edge_id << ", " <<
00067     "No element topology found");
00068 }
00069 
00070 } // unnamed namespace
00071 
00072 Entity & declare_element( BulkData & mesh ,
00073                           Part & part ,
00074                           const EntityId elem_id ,
00075                           const EntityId node_id[] )
00076 {
00077   FEMMetaData & fem_meta = FEMMetaData::get(mesh);
00078   const CellTopologyData * const top = fem_meta.get_cell_topology( part ).getCellTopologyData();
00079 
00080   ThrowErrorMsgIf(top == NULL,
00081                   "Part " << part.name() << " does not have a local topology");
00082 
00083   PartVector empty ;
00084   PartVector add( 1 ); add[0] = & part ;
00085 
00086   const EntityRank entity_rank = fem_meta.element_rank();
00087 
00088   Entity & elem = mesh.declare_entity( entity_rank, elem_id, add );
00089 
00090   const EntityRank node_rank = fem_meta.node_rank();
00091 
00092   for ( unsigned i = 0 ; i < top->node_count ; ++i ) {
00093     //declare node if it doesn't already exist
00094     Entity * node = mesh.get_entity( node_rank , node_id[i]);
00095     if ( NULL == node) {
00096       node = & mesh.declare_entity( node_rank , node_id[i], empty );
00097     }
00098 
00099     mesh.declare_relation( elem , *node , i );
00100   }
00101   return elem ;
00102 }
00103 
00104 Entity & declare_element_side(
00105   Entity & elem ,
00106   Entity & side,
00107   const unsigned local_side_id ,
00108   Part * part )
00109 {
00110   BulkData & mesh = BulkData::get(side);
00111 
00112   verify_declare_element_side(mesh, elem, local_side_id);
00113 
00114   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00115 
00116   ThrowErrorMsgIf( elem_top == NULL,
00117       "Element[" << elem.identifier() << "] has no defined topology" );
00118 
00119   const CellTopologyData * const side_top = elem_top->side[ local_side_id ].topology;
00120 
00121   ThrowErrorMsgIf( side_top == NULL,
00122       "Element[" << elem.identifier() << "], local_side_id = " <<
00123       local_side_id << ", side has no defined topology" );
00124 
00125   const unsigned * const side_node_map = elem_top->side[ local_side_id ].node ;
00126 
00127   PartVector add_parts ;
00128 
00129   if ( part ) { add_parts.push_back( part ); }
00130 
00131   mesh.change_entity_parts(side, add_parts);
00132 
00133   mesh.declare_relation( elem , side , local_side_id );
00134 
00135   PairIterRelation rel = elem.relations( FEMMetaData::NODE_RANK );
00136 
00137   for ( unsigned i = 0 ; i < side_top->node_count ; ++i ) {
00138     Entity & node = * rel[ side_node_map[i] ].entity();
00139     mesh.declare_relation( side , node , i );
00140   }
00141 
00142   return side ;
00143 }
00144 
00145 Entity & declare_element_edge(
00146   Entity & elem ,
00147   Entity & edge,
00148   const unsigned local_edge_id ,
00149   Part * part )
00150 {
00151   BulkData & mesh = BulkData::get(edge);
00152 
00153   //  verify_declare_element_edge(mesh, elem, local_edge_id);
00154 
00155   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00156 
00157   ThrowErrorMsgIf( elem_top == NULL,
00158       "Element[" << elem.identifier() << "] has no defined topology" );
00159 
00160   const CellTopologyData * const edge_top = elem_top->edge[ local_edge_id ].topology;
00161 
00162   ThrowErrorMsgIf( edge_top == NULL,
00163       "Element[" << elem.identifier() << "], local_edge_id = " <<
00164       local_edge_id << ", edge has no defined topology" );
00165 
00166   const unsigned * const edge_node_map = elem_top->edge[ local_edge_id ].node ;
00167 
00168   PartVector add_parts ;
00169 
00170   if ( part ) { add_parts.push_back( part ); }
00171 
00172   mesh.change_entity_parts(edge, add_parts);
00173 
00174   mesh.declare_relation( elem , edge , local_edge_id );
00175 
00176   PairIterRelation rel = elem.relations( FEMMetaData::NODE_RANK );
00177 
00178   for ( unsigned i = 0 ; i < edge_top->node_count ; ++i ) {
00179     Entity & node = * rel[ edge_node_map[i] ].entity();
00180     mesh.declare_relation( edge , node , i );
00181   }
00182 
00183   return edge ;
00184 }
00185 
00186 Entity & declare_element_side(
00187   BulkData & mesh ,
00188   const stk::mesh::EntityId global_side_id ,
00189   Entity & elem ,
00190   const unsigned local_side_id ,
00191   Part * part )
00192 {
00193   verify_declare_element_side(mesh, elem, local_side_id);
00194 
00195   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00196 
00197   ThrowErrorMsgIf( elem_top == NULL,
00198       "Element[" << elem.identifier() << "] has no defined topology");
00199 
00200   const CellTopologyData * const side_top = elem_top->side[ local_side_id ].topology;
00201 
00202   ThrowErrorMsgIf( side_top == NULL,
00203        "Element[" << elem.identifier() << "], local_side_id = " <<
00204        local_side_id << ", side has no defined topology" );
00205 
00206   PartVector empty_parts ;
00207   Entity & side = mesh.declare_entity( side_top->dimension , global_side_id, empty_parts );
00208   return declare_element_side( elem, side, local_side_id, part);
00209 }
00210 
00211 Entity & declare_element_edge(
00212   BulkData & mesh ,
00213   const stk::mesh::EntityId global_edge_id ,
00214   Entity & elem ,
00215   const unsigned local_edge_id ,
00216   Part * part )
00217 {
00218   verify_declare_element_edge(mesh, elem, local_edge_id);
00219 
00220   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00221 
00222   ThrowErrorMsgIf( elem_top == NULL,
00223       "Element[" << elem.identifier() << "] has no defined topology");
00224 
00225 
00226   const CellTopologyData * const edge_top = elem_top->edge[ local_edge_id ].topology;
00227 
00228   ThrowErrorMsgIf( edge_top == NULL,
00229       "Element[" << elem.identifier() << "], local_edge_id = " <<
00230       local_edge_id << ", edge has no defined topology" );
00231 
00232   PartVector empty_parts ;
00233   Entity & edge = mesh.declare_entity( edge_top->dimension , global_edge_id, empty_parts );
00234   return declare_element_edge( elem, edge, local_edge_id, part);
00235 }
00236 
00237 
00238 
00239 const CellTopologyData * get_subcell_nodes(const Entity & entity ,
00240                                            EntityRank subcell_rank ,
00241                                            unsigned subcell_identifier ,
00242                                            EntityVector & subcell_nodes)
00243 {
00244   subcell_nodes.clear();
00245 
00246   // get cell topology
00247   const CellTopologyData* celltopology = get_cell_topology(entity).getCellTopologyData();
00248 
00249   //error checking
00250   {
00251     //no celltopology defined
00252     if (celltopology == NULL) {
00253       return NULL;
00254     }
00255 
00256     // valid ranks fall within the dimension of the cell topology
00257     const bool bad_rank = subcell_rank >= celltopology->dimension;
00258     ThrowInvalidArgMsgIf( bad_rank, "subcell_rank is >= celltopology dimension\n");
00259 
00260     // subcell_identifier must be less than the subcell count
00261     const bool bad_id = subcell_identifier >= celltopology->subcell_count[subcell_rank];
00262     ThrowInvalidArgMsgIf( bad_id,   "subcell_id is >= subcell_count\n");
00263   }
00264 
00265   // Get the cell topology of the subcell
00266   const CellTopologyData * subcell_topology =
00267     celltopology->subcell[subcell_rank][subcell_identifier].topology;
00268 
00269   const int num_nodes_in_subcell = subcell_topology->node_count;
00270 
00271   // For the subcell, get it's local nodes ids
00272   const unsigned* subcell_node_local_ids =
00273     celltopology->subcell[subcell_rank][subcell_identifier].node;
00274 
00275   FEMMetaData & fem_meta = FEMMetaData::get(entity);
00276   const EntityRank node_rank = fem_meta.node_rank();
00277   PairIterRelation node_relations = entity.relations(node_rank);
00278 
00279   subcell_nodes.reserve(num_nodes_in_subcell);
00280 
00281   for (int i = 0; i < num_nodes_in_subcell; ++i ) {
00282     subcell_nodes.push_back( node_relations[subcell_node_local_ids[i]].entity() );
00283   }
00284 
00285   return subcell_topology;
00286 }
00287 
00288 
00289 int get_entity_subcell_id( const Entity & entity ,
00290                            const EntityRank subcell_rank,
00291                            const CellTopologyData * subcell_topology,
00292                            const std::vector<Entity*>& subcell_nodes )
00293 {
00294   const int INVALID_SIDE = -1;
00295 
00296   unsigned num_nodes = subcell_topology->node_count;
00297 
00298   if (num_nodes != subcell_nodes.size()) {
00299     return INVALID_SIDE;
00300   }
00301 
00302   // get topology of elem
00303   const CellTopologyData* entity_topology = get_cell_topology(entity).getCellTopologyData();
00304   if (entity_topology == NULL) {
00305     return INVALID_SIDE;
00306   }
00307 
00308   // get nodal relations for entity
00309   FEMMetaData & fem_meta = FEMMetaData::get(entity);
00310   const EntityRank node_rank = fem_meta.node_rank();
00311   PairIterRelation relations = entity.relations(node_rank);
00312 
00313   const int num_permutations = subcell_topology->permutation_count;
00314 
00315   // Iterate over the subcells of entity...
00316   for (unsigned local_subcell_ordinal = 0;
00317       local_subcell_ordinal < entity_topology->subcell_count[subcell_rank];
00318       ++local_subcell_ordinal) {
00319 
00320     // get topological data for this subcell
00321     const CellTopologyData* curr_subcell_topology =
00322       entity_topology->subcell[subcell_rank][local_subcell_ordinal].topology;
00323 
00324     // If topologies are not the same, there is no way the subcells are the same
00325     if (subcell_topology == curr_subcell_topology) {
00326 
00327       const unsigned* const subcell_node_map = entity_topology->subcell[subcell_rank][local_subcell_ordinal].node;
00328 
00329       // Taking all positive permutations into account, check if this subcell
00330       // has the same nodes as the subcell_nodes argument. Note that this
00331       // implementation preserves the node-order so that we can take
00332       // entity-orientation into account.
00333       for (int p = 0; p < num_permutations; ++p) {
00334 
00335         if (curr_subcell_topology->permutation[p].polarity ==
00336             CELL_PERMUTATION_POLARITY_POSITIVE) {
00337 
00338           const unsigned * const perm_node =
00339             curr_subcell_topology->permutation[p].node ;
00340 
00341           bool all_match = true;
00342           for (unsigned j = 0 ; j < num_nodes; ++j ) {
00343             if (subcell_nodes[j] !=
00344                 relations[subcell_node_map[perm_node[j]]].entity()) {
00345               all_match = false;
00346               break;
00347             }
00348           }
00349 
00350           // all nodes were the same, we have a match
00351           if ( all_match ) {
00352             return local_subcell_ordinal ;
00353           }
00354         }
00355       }
00356     }
00357   }
00358 
00359   return INVALID_SIDE;
00360 }
00361 
00362 bool comm_mesh_counts( BulkData & M ,
00363                        std::vector<size_t> & counts ,
00364                        bool local_flag )
00365 {
00366   const size_t zero = 0 ;
00367 
00368   // Count locally owned entities
00369 
00370   const FEMMetaData & S = FEMMetaData::get(M);
00371   const unsigned entity_rank_count = S.entity_rank_count();
00372   const size_t   comm_count        = entity_rank_count + 1 ;
00373 
00374   std::vector<size_t> local(  comm_count , zero );
00375   std::vector<size_t> global( comm_count , zero );
00376 
00377   ParallelMachine comm = M.parallel();
00378   Part & owns = S.locally_owned_part();
00379 
00380   for ( unsigned i = 0 ; i < entity_rank_count ; ++i ) {
00381     const std::vector<Bucket*> & ks = M.buckets( i );
00382 
00383     std::vector<Bucket*>::const_iterator ik ;
00384 
00385     for ( ik = ks.begin() ; ik != ks.end() ; ++ik ) {
00386       if ( has_superset( **ik , owns ) ) {
00387         local[i] += (*ik)->size();
00388       }
00389     }
00390   }
00391 
00392   local[ entity_rank_count ] = local_flag ;
00393 
00394   stk::all_reduce_sum( comm , & local[0] , & global[0] , comm_count );
00395 
00396   counts.assign( global.begin() , global.begin() + entity_rank_count );
00397 
00398   return 0 < global[ entity_rank_count ] ;
00399 }
00400 
00401 bool element_side_polarity( const Entity & elem ,
00402                             const Entity & side , int local_side_id )
00403 {
00404   // 09/14/10:  TODO:  tscoffe:  Will this work in 1D?
00405   FEMMetaData &fem_meta = FEMMetaData::get(elem);
00406   const bool is_side = side.entity_rank() != fem_meta.edge_rank();
00407   const CellTopologyData * const elem_top = get_cell_topology( elem ).getCellTopologyData();
00408 
00409   const unsigned side_count = ! elem_top ? 0 : (
00410                                 is_side ? elem_top->side_count
00411                                         : elem_top->edge_count );
00412 
00413   ThrowErrorMsgIf( elem_top == NULL,
00414                    "For Element[" << elem.identifier() << "], element has no defined topology");
00415 
00416   ThrowErrorMsgIf( local_side_id < 0 || static_cast<int>(side_count) <= local_side_id,
00417     "For Element[" << elem.identifier() << "], " <<
00418     "side: " << print_entity_key(side) << ", " <<
00419     "local_side_id = " << local_side_id <<
00420     " ; unsupported local_side_id");
00421 
00422   const CellTopologyData * const side_top =
00423     is_side ? elem_top->side[ local_side_id ].topology
00424             : elem_top->edge[ local_side_id ].topology ;
00425 
00426   const unsigned * const side_map =
00427     is_side ? elem_top->side[ local_side_id ].node
00428             : elem_top->edge[ local_side_id ].node ;
00429 
00430   const PairIterRelation elem_nodes = elem.relations( FEMMetaData::NODE_RANK );
00431   const PairIterRelation side_nodes = side.relations( FEMMetaData::NODE_RANK );
00432 
00433   const unsigned n = side_top->node_count;
00434   bool good = false ;
00435   for ( unsigned i = 0 ; !good && i < n ; ++i ) {
00436     good = true;
00437     for ( unsigned j = 0; good && j < n ; ++j ) {
00438       good = side_nodes[(j+i)%n].entity() == elem_nodes[ side_map[j] ].entity();
00439     }
00440   }
00441   return good ;
00442 }
00443 
00444 }
00445 }
00446 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines