MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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