MOAB: Mesh Oriented datABase  (version 5.2.1)
MeshOutputFunctor.cpp
Go to the documentation of this file.
00001 #include "MeshOutputFunctor.hpp"
00002 
00003 #include "SplitVertices.hpp"
00004 #include "moab/ParallelComm.hpp"
00005 #include "RefinerTagManager.hpp"
00006 
00007 #include <iostream>
00008 #include <set>
00009 #include <iterator>
00010 #include <algorithm>
00011 
00012 #undef MB_DEBUG
00013 
00014 namespace moab
00015 {
00016 
00017 MeshOutputFunctor::MeshOutputFunctor( RefinerTagManager* tag_mgr )
00018 {
00019     this->mesh_in         = tag_mgr->get_input_mesh();
00020     this->mesh_out        = tag_mgr->get_output_mesh();
00021     this->input_is_output = ( this->mesh_in == this->mesh_out );
00022     this->tag_manager     = tag_mgr;
00023     this->destination_set = 0;  // don't place output entities in a set by default.
00024 
00025     // When the input mesh and the output mesh are different, this map serves
00026     // as a dictionary from input vertices to output vertices.
00027     this->vertex_map = new SplitVertices< 1 >( this->tag_manager );
00028 
00029     // Hold information about newly-created vertices on subdivided edges and faces.
00030     this->split_vertices.resize( 4 );
00031     this->split_vertices[0] = 0;  // Vertices (0-faces) cannot be split
00032     this->split_vertices[1] = new SplitVertices< 1 >( this->tag_manager );
00033     this->split_vertices[2] = new SplitVertices< 2 >( this->tag_manager );
00034     this->split_vertices[3] = new SplitVertices< 3 >( this->tag_manager );
00035 
00036     // Hold information about newly-created mesh entities (other than split vertices)
00037     // This is necessary in order for global IDs to be assigned consistently across processes.
00038     this->new_entities.resize( 5 );
00039     this->new_entities[0] = new EntitySource( 1, this->tag_manager );
00040     this->new_entities[1] = new EntitySource( 2, this->tag_manager );
00041     this->new_entities[2] = new EntitySource( 3, this->tag_manager );
00042     this->new_entities[3] = new EntitySource( 4, this->tag_manager );
00043     this->new_entities[4] = new EntitySource( 5, this->tag_manager );
00044 }
00045 
00046 MeshOutputFunctor::~MeshOutputFunctor()
00047 {
00048     delete this->vertex_map;
00049     for( int i = 1; i < 4; ++i )
00050         delete this->split_vertices[i];
00051     for( int i = 0; i < 5; ++i )
00052         delete this->new_entities[i];
00053 }
00054 
00055 void MeshOutputFunctor::print_vert_crud( EntityHandle vout, int nvhash, EntityHandle* vhash, const double* vcoords,
00056                                          const void* /*vtags*/ )
00057 {
00058     std::cout << "+ {";
00059     for( int i = 0; i < nvhash; ++i )
00060         std::cout << " " << vhash[i];
00061     std::cout << " } -> " << vout << " ";
00062 
00063     std::cout << "[ " << vcoords[0];
00064     for( int i = 1; i < 6; ++i )
00065         std::cout << ", " << vcoords[i];
00066     std::cout << " ] ";
00067 
00068 #if 0
00069   double* x = (double*)vtags;
00070   int* m = (int*)( (char*)vtags + 2 * sizeof( double ) );
00071   std::cout << "< " << x[0]
00072     << ", " << x[1];
00073   for ( int i = 0; i < 4; ++i )
00074     std::cout << ", " << m[i];
00075 #endif  // 0
00076 
00077     std::cout << " >\n";
00078     // std::cout << "##############################\n";
00079     // this->mesh_out->list_entities( 0, 1 );
00080     // std::cout << "##############################\n";
00081 }
00082 
00083 void MeshOutputFunctor::assign_global_ids( ParallelComm* comm )
00084 {
00085     // First, we must gather the number of entities in each
00086     // partition (for all partitions, not just those resident locally).
00087     int lnparts = this->proc_partition_counts.size();
00088     std::vector< unsigned char > lpdefns;
00089     std::vector< int > lpsizes;
00090     lpdefns.resize( ProcessSet::SHARED_PROC_BYTES * lnparts );
00091     lpsizes.resize( lnparts );
00092 #ifdef MB_DEBUG
00093     std::cout << "**** Partition Counts ****\n";
00094 #endif  // MB_DEBUG
00095     int i = 0;
00096     std::map< ProcessSet, int >::iterator it;
00097     for( it = this->proc_partition_counts.begin(); it != this->proc_partition_counts.end(); ++it, ++i )
00098     {
00099         for( int j = 0; j < ProcessSet::SHARED_PROC_BYTES; ++j )
00100             lpdefns[ProcessSet::SHARED_PROC_BYTES * i + j] = it->first.data()[j];
00101         lpsizes[i] = it->second;
00102 #ifdef MB_DEBUG
00103         std::cout << "Partition " << it->first << ": " << it->second << "\n";
00104 #endif  // MB_DEBUG
00105     }
00106 
00107     if( !comm ) return;
00108 
00109     std::vector< int > nparts;
00110     std::vector< int > dparts;
00111     // unsigned long prank = comm->proc_config().proc_rank();
00112     unsigned long psize = comm->proc_config().proc_size();
00113     nparts.resize( psize );
00114     dparts.resize( psize + 1 );
00115     MPI_Allgather( &lnparts, 1, MPI_INT, &nparts[0], 1, MPI_INT, comm->proc_config().proc_comm() );
00116     // unsigned long ndefs = 0;
00117     for( unsigned long rank = 1; rank <= psize; ++rank )
00118     {
00119         dparts[rank] = nparts[rank - 1] + dparts[rank - 1];
00120 #ifdef MB_DEBUG
00121         std::cout << "Proc " << rank << ": " << nparts[rank - 1] << " partitions, offset: " << dparts[rank] << "\n";
00122 #endif  // MB_DEBUG
00123     }
00124     std::vector< unsigned char > part_defns;
00125     std::vector< int > part_sizes;
00126     part_defns.resize( ProcessSet::SHARED_PROC_BYTES * dparts[psize] );
00127     part_sizes.resize( dparts[psize] );
00128     MPI_Allgatherv( &lpsizes[0], lnparts, MPI_INT, &part_sizes[0], &nparts[0], &dparts[0], MPI_INT,
00129                     comm->proc_config().proc_comm() );
00130     for( unsigned long rank = 0; rank < psize; ++rank )
00131     {
00132         nparts[rank] *= ProcessSet::SHARED_PROC_BYTES;
00133         dparts[rank] *= ProcessSet::SHARED_PROC_BYTES;
00134     }
00135     MPI_Allgatherv( &lpdefns[0], ProcessSet::SHARED_PROC_BYTES * lnparts, MPI_UNSIGNED_CHAR, &part_defns[0], &nparts[0],
00136                     &dparts[0], MPI_UNSIGNED_CHAR, comm->proc_config().proc_comm() );
00137 
00138     // Now that we have the number of new entities in every partition, we
00139     // can deterministically assign the same GID to the same entity even
00140     // when shared across processors because we have an ordering that is
00141     // identical on all processes -- the vertex splits.
00142     for( int j = 0; j < dparts[psize]; ++j )
00143     {
00144         ProcessSet pset( &part_defns[ProcessSet::SHARED_PROC_BYTES * j] );
00145         std::map< ProcessSet, int >::iterator mit = this->proc_partition_counts.find( pset );
00146         if( mit != this->proc_partition_counts.end() )
00147         {
00148 #ifdef MB_DEBUG
00149             std::cout << "Partition " << pset << ( mit->second == part_sizes[j] ? " matches" : " broken" ) << ".\n";
00150 #endif  // MB_DEBUG
00151         }
00152         else
00153         {
00154             this->proc_partition_counts[pset] = part_sizes[j];
00155         }
00156     }
00157 
00158     std::map< ProcessSet, int > gids;
00159     std::map< ProcessSet, int >::iterator pcit;
00160     EntityHandle start_gid = 100;  // FIXME: Get actual maximum GID across all processes and add 1
00161     for( pcit = this->proc_partition_counts.begin(); pcit != this->proc_partition_counts.end(); ++pcit )
00162     {
00163         gids[pcit->first] = start_gid;
00164         start_gid += pcit->second;
00165 #ifdef MB_DEBUG
00166         std::cout << "Partition " << pcit->first << ": " << pcit->second << " # [" << gids[pcit->first] << "]\n";
00167 #endif  // MB_DEBUG
00168     }
00169     std::vector< SplitVerticesBase* >::iterator vit;
00170     vit = this->split_vertices.begin();
00171     ++vit;  // Skip split_vertices[0] since it's empty.
00172     ++vit;  // Skip split_vertices[1] since those entries already have global IDs... they exist in
00173             // the input mesh.
00174     for( /* skip */; vit != this->split_vertices.end(); ++vit )
00175     {
00176         ( *vit )->assign_global_ids( gids );
00177     }
00178     std::vector< EntitySource* >::iterator sit;
00179     for( sit = this->new_entities.begin(); sit != this->new_entities.end(); ++sit )
00180     {
00181         if( *sit ) ( *sit )->assign_global_ids( gids );
00182     }
00183 }
00184 
00185 void MeshOutputFunctor::exchange_handles( ParallelComm* ) {}
00186 
00187 void MeshOutputFunctor::assign_tags( EntityHandle vhandle, const void* vtags )
00188 {
00189     if( !vhandle ) return;  // Ignore bad vertices
00190 
00191     int num_tags = this->tag_manager->get_number_of_vertex_tags();
00192     Tag tag_handle;
00193     int tag_offset;
00194     for( int i = 0; i < num_tags; ++i )
00195     {
00196         this->tag_manager->get_output_vertex_tag( i, tag_handle, tag_offset );
00197         this->mesh_out->tag_set_data( tag_handle, &vhandle, 1, vtags );
00198     }
00199 }
00200 
00201 EntityHandle MeshOutputFunctor::map_vertex( EntityHandle vhash, const double* vcoords, const void* vtags )
00202 {
00203     if( this->input_is_output )
00204     {  // Don't copy the original vertex!
00205 #ifdef MB_DEBUG
00206         this->print_vert_crud( vhash, 1, &vhash, vcoords, vtags );
00207 #endif  // MB_DEBUG
00208         return vhash;
00209     }
00210     EntityHandle vertex_handle;
00211     bool newly_created =
00212         this->vertex_map->find_or_create( &vhash, vcoords, vertex_handle, this->proc_partition_counts, false );
00213     if( newly_created )
00214     {
00215         std::vector< int > gid;
00216         this->assign_tags( vertex_handle, vtags );
00217         if( this->tag_manager->get_input_gids( 1, &vhash, gid ) == MB_SUCCESS )
00218         { this->tag_manager->set_gid( vertex_handle, gid[0] ); }
00219     }
00220     if( !vertex_handle ) { std::cerr << "Could not insert vertex into new mesh!\n"; }
00221 #ifdef MB_DEBUG
00222     this->print_vert_crud( vertex_handle, 1, &vhash, vcoords, vtags );
00223     std::cout << "\nMap vert: " << vhash << " to: " << vertex_handle << "\n";
00224 #endif  // MB_DEBUG
00225     return vertex_handle;
00226 }
00227 
00228 EntityHandle MeshOutputFunctor::operator()( int nvhash, EntityHandle* vhash, const double* vcoords, const void* vtags )
00229 {
00230     EntityHandle vertex_handle;
00231     if( nvhash < 4 )
00232     {
00233         bool newly_created = this->split_vertices[nvhash]->find_or_create( vhash, vcoords, vertex_handle,
00234                                                                            this->proc_partition_counts, true );
00235         if( newly_created ) { this->assign_tags( vertex_handle, vtags ); }
00236         if( !vertex_handle ) { std::cerr << "Could not insert mid-edge vertex!\n"; }
00237 #ifdef MB_DEBUG
00238         std::cout << "(-" << nvhash << "-) ";
00239         this->print_vert_crud( vertex_handle, nvhash, vhash, vcoords, vtags );
00240 #endif  // MB_DEBUG
00241     }
00242     else
00243     {
00244         vertex_handle = 0;
00245         std::cerr << "Not handling splits on faces with " << nvhash << " corners yet.\n";
00246     }
00247     return vertex_handle;
00248 }
00249 
00250 void MeshOutputFunctor::operator()( EntityHandle h )
00251 {
00252 #ifdef MB_DEBUG
00253     std::cout << h << " ";
00254 #endif  // MB_DEBUG
00255     if( !this->input_is_output )
00256     {
00257         // FIXME: Copy to output mesh
00258     }
00259     this->elem_vert.push_back( h );
00260 }
00261 
00262 void MeshOutputFunctor::operator()( EntityType etyp )
00263 {
00264     EntityHandle elem_handle;
00265     int nconn          = this->elem_vert.size();
00266     bool newly_created = this->new_entities[nconn]->create_element( etyp, nconn, &this->elem_vert[0], elem_handle,
00267                                                                     this->proc_partition_counts );
00268     if( newly_created )
00269     {
00270 #ifdef MB_DEBUG
00271         std::cout << " *** ";
00272 #endif  // MB_DEBUG
00273         // FIXME: Handle tag assignment for elements as well as vertices
00274         this->tag_manager->assign_element_tags( elem_handle );
00275     }
00276 #ifdef MB_DEBUG
00277     std::cout << "---------> " << elem_handle << " ( " << etyp << " )\n\n";
00278 #endif  // MB_DEBUG
00279     this->elem_vert.clear();
00280 }
00281 
00282 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines