MOAB: Mesh Oriented datABase  (version 5.4.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,
00056                                          int nvhash,
00057                                          EntityHandle* vhash,
00058                                          const double* vcoords,
00059                                          const void* /*vtags*/ )
00060 {
00061     std::cout << "+ {";
00062     for( int i = 0; i < nvhash; ++i )
00063         std::cout << " " << vhash[i];
00064     std::cout << " } -> " << vout << " ";
00065 
00066     std::cout << "[ " << vcoords[0];
00067     for( int i = 1; i < 6; ++i )
00068         std::cout << ", " << vcoords[i];
00069     std::cout << " ] ";
00070 
00071 #if 0
00072   double* x = (double*)vtags;
00073   int* m = (int*)( (char*)vtags + 2 * sizeof( double ) );
00074   std::cout << "< " << x[0]
00075     << ", " << x[1];
00076   for ( int i = 0; i < 4; ++i )
00077     std::cout << ", " << m[i];
00078 #endif  // 0
00079 
00080     std::cout << " >\n";
00081     // std::cout << "##############################\n";
00082     // this->mesh_out->list_entities( 0, 1 );
00083     // std::cout << "##############################\n";
00084 }
00085 
00086 void MeshOutputFunctor::assign_global_ids( ParallelComm* comm )
00087 {
00088     // First, we must gather the number of entities in each
00089     // partition (for all partitions, not just those resident locally).
00090     int lnparts = this->proc_partition_counts.size();
00091     std::vector< unsigned char > lpdefns;
00092     std::vector< int > lpsizes;
00093     lpdefns.resize( ProcessSet::SHARED_PROC_BYTES * lnparts );
00094     lpsizes.resize( lnparts );
00095 #ifdef MB_DEBUG
00096     std::cout << "**** Partition Counts ****\n";
00097 #endif  // MB_DEBUG
00098     int i = 0;
00099     std::map< ProcessSet, int >::iterator it;
00100     for( it = this->proc_partition_counts.begin(); it != this->proc_partition_counts.end(); ++it, ++i )
00101     {
00102         for( int j = 0; j < ProcessSet::SHARED_PROC_BYTES; ++j )
00103             lpdefns[ProcessSet::SHARED_PROC_BYTES * i + j] = it->first.data()[j];
00104         lpsizes[i] = it->second;
00105 #ifdef MB_DEBUG
00106         std::cout << "Partition " << it->first << ": " << it->second << "\n";
00107 #endif  // MB_DEBUG
00108     }
00109 
00110     if( !comm ) return;
00111 
00112     std::vector< int > nparts;
00113     std::vector< int > dparts;
00114     // unsigned long prank = comm->proc_config().proc_rank();
00115     unsigned long psize = comm->proc_config().proc_size();
00116     nparts.resize( psize );
00117     dparts.resize( psize + 1 );
00118     MPI_Allgather( &lnparts, 1, MPI_INT, &nparts[0], 1, MPI_INT, comm->proc_config().proc_comm() );
00119     // unsigned long ndefs = 0;
00120     for( unsigned long rank = 1; rank <= psize; ++rank )
00121     {
00122         dparts[rank] = nparts[rank - 1] + dparts[rank - 1];
00123 #ifdef MB_DEBUG
00124         std::cout << "Proc " << rank << ": " << nparts[rank - 1] << " partitions, offset: " << dparts[rank] << "\n";
00125 #endif  // MB_DEBUG
00126     }
00127     std::vector< unsigned char > part_defns;
00128     std::vector< int > part_sizes;
00129     part_defns.resize( ProcessSet::SHARED_PROC_BYTES * dparts[psize] );
00130     part_sizes.resize( dparts[psize] );
00131     MPI_Allgatherv( &lpsizes[0], lnparts, MPI_INT, &part_sizes[0], &nparts[0], &dparts[0], MPI_INT,
00132                     comm->proc_config().proc_comm() );
00133     for( unsigned long rank = 0; rank < psize; ++rank )
00134     {
00135         nparts[rank] *= ProcessSet::SHARED_PROC_BYTES;
00136         dparts[rank] *= ProcessSet::SHARED_PROC_BYTES;
00137     }
00138     MPI_Allgatherv( &lpdefns[0], ProcessSet::SHARED_PROC_BYTES * lnparts, MPI_UNSIGNED_CHAR, &part_defns[0], &nparts[0],
00139                     &dparts[0], MPI_UNSIGNED_CHAR, comm->proc_config().proc_comm() );
00140 
00141     // Now that we have the number of new entities in every partition, we
00142     // can deterministically assign the same GID to the same entity even
00143     // when shared across processors because we have an ordering that is
00144     // identical on all processes -- the vertex splits.
00145     for( int j = 0; j < dparts[psize]; ++j )
00146     {
00147         ProcessSet pset( &part_defns[ProcessSet::SHARED_PROC_BYTES * j] );
00148         std::map< ProcessSet, int >::iterator mit = this->proc_partition_counts.find( pset );
00149         if( mit != this->proc_partition_counts.end() )
00150         {
00151 #ifdef MB_DEBUG
00152             std::cout << "Partition " << pset << ( mit->second == part_sizes[j] ? " matches" : " broken" ) << ".\n";
00153 #endif  // MB_DEBUG
00154         }
00155         else
00156         {
00157             this->proc_partition_counts[pset] = part_sizes[j];
00158         }
00159     }
00160 
00161     std::map< ProcessSet, int > gids;
00162     std::map< ProcessSet, int >::iterator pcit;
00163     EntityHandle start_gid = 100;  // FIXME: Get actual maximum GID across all processes and add 1
00164     for( pcit = this->proc_partition_counts.begin(); pcit != this->proc_partition_counts.end(); ++pcit )
00165     {
00166         gids[pcit->first] = start_gid;
00167         start_gid += pcit->second;
00168 #ifdef MB_DEBUG
00169         std::cout << "Partition " << pcit->first << ": " << pcit->second << " # [" << gids[pcit->first] << "]\n";
00170 #endif  // MB_DEBUG
00171     }
00172     std::vector< SplitVerticesBase* >::iterator vit;
00173     vit = this->split_vertices.begin();
00174     ++vit;  // Skip split_vertices[0] since it's empty.
00175     ++vit;  // Skip split_vertices[1] since those entries already have global IDs... they exist in
00176             // the input mesh.
00177     for( /* skip */; vit != this->split_vertices.end(); ++vit )
00178     {
00179         ( *vit )->assign_global_ids( gids );
00180     }
00181     std::vector< EntitySource* >::iterator sit;
00182     for( sit = this->new_entities.begin(); sit != this->new_entities.end(); ++sit )
00183     {
00184         if( *sit ) ( *sit )->assign_global_ids( gids );
00185     }
00186 }
00187 
00188 void MeshOutputFunctor::exchange_handles( ParallelComm* ) {}
00189 
00190 void MeshOutputFunctor::assign_tags( EntityHandle vhandle, const void* vtags )
00191 {
00192     if( !vhandle ) return;  // Ignore bad vertices
00193 
00194     int num_tags = this->tag_manager->get_number_of_vertex_tags();
00195     Tag tag_handle;
00196     int tag_offset;
00197     for( int i = 0; i < num_tags; ++i )
00198     {
00199         this->tag_manager->get_output_vertex_tag( i, tag_handle, tag_offset );
00200         this->mesh_out->tag_set_data( tag_handle, &vhandle, 1, vtags );
00201     }
00202 }
00203 
00204 EntityHandle MeshOutputFunctor::map_vertex( EntityHandle vhash, const double* vcoords, const void* vtags )
00205 {
00206     if( this->input_is_output )
00207     {  // Don't copy the original vertex!
00208 #ifdef MB_DEBUG
00209         this->print_vert_crud( vhash, 1, &vhash, vcoords, vtags );
00210 #endif  // MB_DEBUG
00211         return vhash;
00212     }
00213     EntityHandle vertex_handle;
00214     bool newly_created =
00215         this->vertex_map->find_or_create( &vhash, vcoords, vertex_handle, this->proc_partition_counts, false );
00216     if( newly_created )
00217     {
00218         std::vector< int > gid;
00219         this->assign_tags( vertex_handle, vtags );
00220         if( this->tag_manager->get_input_gids( 1, &vhash, gid ) == MB_SUCCESS )
00221         {
00222             this->tag_manager->set_gid( vertex_handle, gid[0] );
00223         }
00224     }
00225     if( !vertex_handle )
00226     {
00227         std::cerr << "Could not insert vertex into new mesh!\n";
00228     }
00229 #ifdef MB_DEBUG
00230     this->print_vert_crud( vertex_handle, 1, &vhash, vcoords, vtags );
00231     std::cout << "\nMap vert: " << vhash << " to: " << vertex_handle << "\n";
00232 #endif  // MB_DEBUG
00233     return vertex_handle;
00234 }
00235 
00236 EntityHandle MeshOutputFunctor::operator()( int nvhash, EntityHandle* vhash, const double* vcoords, const void* vtags )
00237 {
00238     EntityHandle vertex_handle;
00239     if( nvhash < 4 )
00240     {
00241         bool newly_created = this->split_vertices[nvhash]->find_or_create( vhash, vcoords, vertex_handle,
00242                                                                            this->proc_partition_counts, true );
00243         if( newly_created )
00244         {
00245             this->assign_tags( vertex_handle, vtags );
00246         }
00247         if( !vertex_handle )
00248         {
00249             std::cerr << "Could not insert mid-edge vertex!\n";
00250         }
00251 #ifdef MB_DEBUG
00252         std::cout << "(-" << nvhash << "-) ";
00253         this->print_vert_crud( vertex_handle, nvhash, vhash, vcoords, vtags );
00254 #endif  // MB_DEBUG
00255     }
00256     else
00257     {
00258         vertex_handle = 0;
00259         std::cerr << "Not handling splits on faces with " << nvhash << " corners yet.\n";
00260     }
00261     return vertex_handle;
00262 }
00263 
00264 void MeshOutputFunctor::operator()( EntityHandle h )
00265 {
00266 #ifdef MB_DEBUG
00267     std::cout << h << " ";
00268 #endif  // MB_DEBUG
00269     if( !this->input_is_output )
00270     {
00271         // FIXME: Copy to output mesh
00272     }
00273     this->elem_vert.push_back( h );
00274 }
00275 
00276 void MeshOutputFunctor::operator()( EntityType etyp )
00277 {
00278     EntityHandle elem_handle;
00279     int nconn          = this->elem_vert.size();
00280     bool newly_created = this->new_entities[nconn]->create_element( etyp, nconn, &this->elem_vert[0], elem_handle,
00281                                                                     this->proc_partition_counts );
00282     if( newly_created )
00283     {
00284 #ifdef MB_DEBUG
00285         std::cout << " *** ";
00286 #endif  // MB_DEBUG 
00287         // FIXME: Handle tag assignment for elements as well as vertices
00288         this->tag_manager->assign_element_tags( elem_handle );
00289     }
00290 #ifdef MB_DEBUG
00291     std::cout << "---------> " << elem_handle << " ( " << etyp << " )\n\n";
00292 #endif  // MB_DEBUG
00293     this->elem_vert.clear();
00294 }
00295 
00296 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines