MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* 00002 * Intx2Mesh.cpp 00003 * 00004 * Created on: Oct 2, 2012 00005 */ 00006 00007 #include "moab/IntxMesh/Intx2Mesh.hpp" 00008 #ifdef MOAB_HAVE_MPI 00009 #include "moab/ParallelComm.hpp" 00010 #include "MBParallelConventions.h" 00011 #include "moab/ParallelMergeMesh.hpp" 00012 #endif /* MOAB_HAVE_MPI */ 00013 #include "MBTagConventions.hpp" 00014 // this is for DBL_MAX 00015 #include <float.h> 00016 #include <queue> 00017 #include <sstream> 00018 #include "moab/GeomUtil.hpp" 00019 #include "moab/AdaptiveKDTree.hpp" 00020 00021 namespace moab 00022 { 00023 00024 #ifdef ENABLE_DEBUG 00025 int Intx2Mesh::dbg_1 = 0; 00026 #endif 00027 00028 Intx2Mesh::Intx2Mesh( Interface* mbimpl ) 00029 : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ), 00030 countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), tgtConn( NULL ), 00031 srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ), my_rank( 0 ) 00032 #ifdef MOAB_HAVE_MPI 00033 , 00034 parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL ) 00035 #endif 00036 , 00037 max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 ) 00038 { 00039 gid = mbimpl->globalId_tag(); 00040 } 00041 00042 Intx2Mesh::~Intx2Mesh() 00043 { 00044 // TODO Auto-generated destructor stub 00045 #ifdef MOAB_HAVE_MPI 00046 if( remote_cells ) 00047 { 00048 delete remote_cells; 00049 remote_cells = NULL; 00050 } 00051 #endif 00052 } 00053 ErrorCode Intx2Mesh::FindMaxEdgesInSet( EntityHandle eset, int& max_edges ) 00054 { 00055 Range cells; 00056 ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval ); 00057 00058 max_edges = 0; // can be 0 for point clouds 00059 for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ ) 00060 { 00061 EntityHandle cell = *cit; 00062 const EntityHandle* conn4; 00063 int nnodes = 3; 00064 rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" ); 00065 if( nnodes > max_edges ) max_edges = nnodes; 00066 } 00067 // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be 00068 // global) but it is better to be consistent 00069 #ifdef MOAB_HAVE_MPI 00070 if( parcomm ) 00071 { 00072 int local_max_edges = max_edges; 00073 // now reduce max_edges over all processors 00074 int mpi_err = 00075 MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() ); 00076 if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 00077 } 00078 #endif 00079 00080 return MB_SUCCESS; 00081 } 00082 ErrorCode Intx2Mesh::FindMaxEdges( EntityHandle set1, EntityHandle set2 ) 00083 { 00084 ErrorCode rval = FindMaxEdgesInSet( set1, max_edges_1 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 1" ); 00085 rval = FindMaxEdgesInSet( set2, max_edges_2 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 2" ); 00086 00087 return MB_SUCCESS; 00088 } 00089 00090 ErrorCode Intx2Mesh::createTags() 00091 { 00092 if( tgtParentTag ) mb->tag_delete( tgtParentTag ); 00093 if( srcParentTag ) mb->tag_delete( srcParentTag ); 00094 if( countTag ) mb->tag_delete( countTag ); 00095 00096 unsigned char def_data_bit = 0; // unused by default 00097 // maybe the tgt tag is better to be deleted every time, and recreated; 00098 // or is it easy to set all values to something again? like 0? 00099 ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "can't get tgt flag tag" ); 00100 00101 // create tgt edges if they do not exist yet; so when they are looked upon, they are found 00102 // this is the only call that is potentially NlogN, in the whole method 00103 rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" ); 00104 00105 // now, create a map from each edge to a list of potential new nodes on a tgt edge 00106 // this memory has to be cleaned up 00107 // change it to a vector, and use the index in range of tgt edges 00108 int indx = 0; 00109 extraNodesVec.resize( TgtEdges.size() ); 00110 for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ ) 00111 { 00112 std::vector< EntityHandle >* nv = new std::vector< EntityHandle >; 00113 extraNodesVec[indx] = nv; 00114 } 00115 00116 int defaultInt = -1; 00117 00118 rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT, 00119 &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" ); 00120 00121 rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT, 00122 &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" ); 00123 00124 rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" ); 00125 00126 // for each cell in set 1, determine its neigh in set 1 (could be null too) 00127 // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0) 00128 rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" ); 00129 rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 2" ); 00130 00131 // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for 00132 // them each time edges were for sure created before (tgtEdges) 00133 std::vector< EntityHandle > zeroh( max_edges_2, 0 ); 00134 // if we have a tag with this name, it could be of a different size, so delete it 00135 rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag ); 00136 if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag ); 00137 rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag, 00138 MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" ); 00139 for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ ) 00140 { 00141 EntityHandle tgtCell = *rit; 00142 int num_nodes = 0; 00143 rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" ); 00144 // account for padded polygons 00145 while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 ) 00146 num_nodes--; 00147 00148 int i = 0; 00149 for( i = 0; i < num_nodes; i++ ) 00150 { 00151 EntityHandle v[2] = { tgtConn[i], 00152 tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons 00153 std::vector< EntityHandle > adj_entities; 00154 rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT ); 00155 if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error 00156 zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes 00157 // also, even if number of edges is less than max_edges_2, they will be ignored, even if 00158 // the tag is dense 00159 } 00160 // now set the value of the tag 00161 rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" ); 00162 } 00163 return MB_SUCCESS; 00164 } 00165 00166 ErrorCode Intx2Mesh::DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag ) 00167 { 00168 Range cells; 00169 ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells in set" ); 00170 00171 std::vector< EntityHandle > neighbors( max_edges ); 00172 std::vector< EntityHandle > zeroh( max_edges, 0 ); 00173 // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh, 00174 // one on src mesh 00175 rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create neighbors tag" ); 00176 00177 for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ ) 00178 { 00179 EntityHandle cell = *cit; 00180 int nnodes = 3; 00181 // will get the nnodes ordered neighbors; 00182 // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1, 00183 const EntityHandle* conn4; 00184 rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" ); 00185 int nsides = nnodes; 00186 // account for possible padded polygons 00187 while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 ) 00188 nsides--; 00189 00190 for( int i = 0; i < nsides; i++ ) 00191 { 00192 EntityHandle v[2]; 00193 v[0] = conn4[i]; 00194 v[1] = conn4[( i + 1 ) % nsides]; 00195 // get all cells adjacent to these 2 vertices on the edge 00196 std::vector< EntityHandle > adjcells; 00197 std::vector< EntityHandle > cellsInSet; 00198 rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't adjacency to 2 verts" ); 00199 // now look for the cells contained in the input set; 00200 // the input set should be a correct mesh, not overlapping cells, and manifold 00201 size_t siz = adjcells.size(); 00202 for( size_t j = 0; j < siz; j++ ) 00203 if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] ); 00204 siz = cellsInSet.size(); 00205 00206 if( siz > 2 ) 00207 { 00208 std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() ) 00209 << "\n";MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" ); // non-manifold 00210 } 00211 if( siz == 1 ) 00212 { 00213 // it must be the border of the input mesh; 00214 neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border 00215 // borders do not appear for a sphere in serial, but they do appear for 00216 // parallel processing anyway 00217 continue; 00218 } 00219 // here siz ==2, it is either the first or second 00220 if( cell == cellsInSet[0] ) 00221 neighbors[i] = cellsInSet[1]; 00222 else 00223 neighbors[i] = cellsInSet[0]; 00224 } 00225 // fill the rest with 0 00226 for( int i = nsides; i < max_edges; i++ ) 00227 neighbors[i] = 0; 00228 // now simply set the neighbors tag; the last few positions will not be used, but for 00229 // simplicity will keep them all (MAXEDGES) 00230 rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" ); 00231 } 00232 return MB_SUCCESS; 00233 } 00234 00235 // slow interface; this will not do the advancing front trick 00236 // some are triangles, some are quads, some are polygons ... 00237 ErrorCode Intx2Mesh::intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet ) 00238 { 00239 ErrorCode rval; 00240 mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc 00241 mbs2 = mbset2; 00242 outSet = outputSet; 00243 rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval ); 00244 rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval ); 00245 // from create tags, copy relevant ones 00246 if( tgtParentTag ) mb->tag_delete( tgtParentTag ); 00247 if( srcParentTag ) mb->tag_delete( srcParentTag ); 00248 if( countTag ) mb->tag_delete( countTag ); 00249 00250 // create tgt edges if they do not exist yet; so when they are looked upon, they are found 00251 // this is the only call that is potentially NlogN, in the whole method 00252 rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" ); 00253 00254 int indx = 0; 00255 extraNodesVec.resize( TgtEdges.size() ); 00256 for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ ) 00257 { 00258 std::vector< EntityHandle >* nv = new std::vector< EntityHandle >; 00259 extraNodesVec[indx] = nv; 00260 } 00261 00262 int defaultInt = -1; 00263 00264 rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT, 00265 &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" ); 00266 00267 rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT, 00268 &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" ); 00269 00270 rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" ); 00271 00272 // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for 00273 // them each time edges were for sure created before (tgtEdges) 00274 // if we have a tag with this name, it could be of a different size, so delete it 00275 rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag ); 00276 if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag ); 00277 std::vector< EntityHandle > zeroh( max_edges_2, 0 ); 00278 rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag, 00279 MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" ); 00280 for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ ) 00281 { 00282 EntityHandle tgtCell = *rit; 00283 int num_nodes = 0; 00284 rval = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get tgt conn" ); 00285 // account for padded polygons 00286 while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 ) 00287 num_nodes--; 00288 00289 int i = 0; 00290 for( i = 0; i < num_nodes; i++ ) 00291 { 00292 EntityHandle v[2] = { tgtConn[i], 00293 tgtConn[( i + 1 ) % num_nodes] }; // this is fine even for padded polygons 00294 std::vector< EntityHandle > adj_entities; 00295 rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT ); 00296 if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval; // get out , big error 00297 zeroh[i] = adj_entities[0]; // should be only one edge between 2 nodes 00298 // also, even if number of edges is less than max_edges_2, they will be ignored, even if 00299 // the tag is dense 00300 } 00301 // now set the value of the tag 00302 rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" ); 00303 } 00304 00305 // create the kd tree on source cells, and intersect all targets in an expensive loop 00306 // build a kd tree with the rs1 (source) cells 00307 AdaptiveKDTree kd( mb ); 00308 EntityHandle tree_root = 0; 00309 rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval ); 00310 00311 // find out max edge on source mesh; 00312 double max_length = 0; 00313 { 00314 std::vector< double > coords; 00315 coords.resize( 3 * max_edges_1 ); 00316 for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ ) 00317 { 00318 const EntityHandle* conn = NULL; 00319 int nnodes; 00320 rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" ); 00321 while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 ) 00322 nnodes--; 00323 rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coordinates" ); 00324 for( int j = 0; j < nnodes; j++ ) 00325 { 00326 int next = ( j + 1 ) % nnodes; 00327 double leng; 00328 leng = ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) + 00329 ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) + 00330 ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] ); 00331 leng = sqrt( leng ); 00332 if( leng > max_length ) max_length = leng; 00333 } 00334 } 00335 } 00336 // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :( 00337 double tolerance = 1.e-15; 00338 if( max_length < 1. ) 00339 { 00340 // basically, the sag for an arc of length max_length on a circle of radius 1 00341 tolerance = 1. - sqrt( 1 - max_length * max_length / 4 ); 00342 tolerance = 3 * tolerance; // why ? 00343 if( !my_rank ) 00344 std::cout << " max edge length: " << max_length << " tolerance for kd tree: " << tolerance << "\n"; 00345 } 00346 for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it ) 00347 { 00348 EntityHandle tcell = *it; 00349 // find vertex positions 00350 const EntityHandle* conn = NULL; 00351 int nnodes = 0; 00352 rval = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval ); 00353 // find leaves close to those positions 00354 double areaTgtCell = setup_tgt_cell( tcell, nnodes ); // this is the area in the gnomonic plane 00355 double recoveredArea = 0; 00356 std::vector< double > positions; 00357 positions.resize( nnodes * 3 ); 00358 rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval ); 00359 00360 // distance to search will be based on average edge length 00361 double av_len = 0; 00362 for( int k = 0; k < nnodes; k++ ) 00363 { 00364 int ik = ( k + 1 ) % nnodes; 00365 double len1 = 0; 00366 for( int j = 0; j < 3; j++ ) 00367 { 00368 double len2 = positions[3 * k + j] - positions[3 * ik + j]; 00369 len1 += len2 * len2; 00370 } 00371 av_len += sqrt( len1 ); 00372 } 00373 if( nnodes > 0 ) av_len /= nnodes; 00374 // find leaves within a distance from each vertex of target 00375 // in those leaves, collect all cells; we will try for an intx in there 00376 Range close_source_cells; 00377 std::vector< EntityHandle > leaves; 00378 for( int i = 0; i < nnodes; i++ ) 00379 { 00380 00381 leaves.clear(); 00382 rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval ); 00383 00384 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j ) 00385 { 00386 Range tmp; 00387 rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval ); 00388 00389 close_source_cells.merge( tmp.begin(), tmp.end() ); 00390 } 00391 } 00392 #ifdef VERBOSE 00393 if( close_source_cells.empty() ) 00394 { 00395 std::cout << " there are no close source cells to target cell " << tcell << " id from handle " 00396 << mb->id_from_handle( tcell ) << "\n"; 00397 } 00398 #endif 00399 for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 ) 00400 { 00401 EntityHandle startSrc = *it2; 00402 double area = 0; 00403 // if area is > 0 , we have intersections 00404 double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon 00405 // 00406 int nP = 0; 00407 int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first 00408 int nsTgt, nsSrc; 00409 rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval ); 00410 if( area > 0 ) 00411 { 00412 if( nP > 1 ) 00413 { // this will also construct triangles/polygons in the new mesh, if needed 00414 rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval ); 00415 } 00416 recoveredArea += area; 00417 } 00418 } 00419 recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fract 00420 } 00421 // before cleaning up , we need to settle the position of the intersection points 00422 // on the boundary edges 00423 // this needs to be collective, so we should maybe wait something 00424 #ifdef MOAB_HAVE_MPI 00425 rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" ); 00426 #endif 00427 00428 this->clean(); 00429 return MB_SUCCESS; 00430 } 00431 // main interface; this will do the advancing front trick 00432 // some are triangles, some are quads, some are polygons ... 00433 ErrorCode Intx2Mesh::intersect_meshes( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet ) 00434 { 00435 ErrorCode rval; 00436 mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc 00437 mbs2 = mbset2; 00438 outSet = outputSet; 00439 #ifdef VERBOSE 00440 std::stringstream ffs, fft; 00441 ffs << "source_rank0" << my_rank << ".vtk"; 00442 rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval ); 00443 fft << "target_rank0" << my_rank << ".vtk"; 00444 rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval ); 00445 00446 #endif 00447 // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler) 00448 00449 EntityHandle startSrc = 0, startTgt = 0; 00450 00451 rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval ); 00452 rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval ); 00453 // std::cout << "rs1.size() = " << rs1.size() << " and rs2.size() = " << rs2.size() << "\n"; 00454 // std::cout.flush(); 00455 00456 createTags(); // will also determine max_edges_1, max_edges_2 (for src and tgt meshes) 00457 00458 Range rs22 = rs2; // a copy of the initial range; we will remove from it elements as we 00459 // advance ; rs2 is needed for marking the polygon to the tgt parent 00460 00461 // create the local kdd tree with source elements; will use it to search 00462 // more efficiently for the seeds in advancing front; 00463 // some of the target cells will not be covered by source cells, and they need to be eliminated 00464 // early from contention 00465 00466 // build a kd tree with the rs1 (source) cells 00467 AdaptiveKDTree kd( mb ); 00468 EntityHandle tree_root = 0; 00469 rval = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval ); 00470 00471 while( !rs22.empty() ) 00472 { 00473 #if defined( ENABLE_DEBUG ) || defined( VERBOSE ) 00474 if( rs22.size() < rs2.size() ) 00475 { 00476 std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting 00477 << "\n"; 00478 std::stringstream ffo; 00479 ffo << "file0" << counting << "rank0" << my_rank << ".vtk"; 00480 rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval ); 00481 } 00482 #endif 00483 bool seedFound = false; 00484 for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it ) 00485 { 00486 startTgt = *it; 00487 int found = 0; 00488 // find vertex positions 00489 const EntityHandle* conn = NULL; 00490 int nnodes = 0; 00491 rval = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval ); 00492 // find leaves close to those positions 00493 std::vector< double > positions; 00494 positions.resize( nnodes * 3 ); 00495 rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval ); 00496 // find leaves within a distance from each vertex of target 00497 // in those leaves, collect all cells; we will try for an intx in there, instead of 00498 // looping over all rs1 cells, as before 00499 Range close_source_cells; 00500 std::vector< EntityHandle > leaves; 00501 for( int i = 0; i < nnodes; i++ ) 00502 { 00503 00504 leaves.clear(); 00505 rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval ); 00506 00507 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j ) 00508 { 00509 Range tmp; 00510 rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval ); 00511 00512 close_source_cells.merge( tmp.begin(), tmp.end() ); 00513 } 00514 } 00515 00516 for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 ) 00517 { 00518 startSrc = *it2; 00519 double area = 0; 00520 // if area is > 0 , we have intersections 00521 double P[10 * MAXEDGES]; // max 8 intx points + 8 more in the polygon 00522 // 00523 int nP = 0; 00524 int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first 00525 int nsTgt, nsSrc; 00526 rval = 00527 computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval ); 00528 if( area > 0 ) 00529 { 00530 found = 1; 00531 seedFound = true; 00532 break; // found 2 elements that intersect; these will be the seeds 00533 } 00534 } 00535 if( found ) 00536 break; 00537 else 00538 { 00539 #if defined( VERBOSE ) 00540 std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt ) 00541 << " not intx with any source\n"; 00542 #endif 00543 rs22.erase( startTgt ); 00544 } 00545 } 00546 if( !seedFound ) continue; // continue while(!rs22.empty()) 00547 00548 std::queue< EntityHandle > srcQueue; // these are corresponding to Ta, 00549 srcQueue.push( startSrc ); 00550 std::queue< EntityHandle > tgtQueue; 00551 tgtQueue.push( startTgt ); 00552 00553 Range toResetSrcs; // will be used to reset src flags for every tgt quad 00554 // processed 00555 00556 /*if (my_rank==0) 00557 dbg_1 = 1;*/ 00558 unsigned char used = 1; 00559 // mark the start tgt quad as used, so it will not come back again 00560 rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval ); 00561 while( !tgtQueue.empty() ) 00562 { 00563 // flags for the side : 0 means a src cell not found on side 00564 // a paired src not found yet for the neighbors of tgt 00565 Range nextSrc[MAXEDGES]; // there are new ranges of possible next src cells for 00566 // seeding the side j of tgt cell 00567 00568 EntityHandle currentTgt = tgtQueue.front(); 00569 tgtQueue.pop(); 00570 int nsidesTgt; // will be initialized now 00571 double areaTgtCell = setup_tgt_cell( currentTgt, nsidesTgt ); // this is the area in the gnomonic plane 00572 double recoveredArea = 0; 00573 // get the neighbors of tgt, and if they are solved already, do not bother with that 00574 // side of tgt 00575 EntityHandle tgtNeighbors[MAXEDGES] = { 0 }; 00576 rval = mb->tag_get_data( tgtNeighTag, ¤tTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" ); 00577 #ifdef ENABLE_DEBUG 00578 if( dbg_1 ) 00579 { 00580 std::cout << "Next: neighbors for current tgt "; 00581 for( int kk = 0; kk < nsidesTgt; kk++ ) 00582 { 00583 if( tgtNeighbors[kk] > 0 ) 00584 std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " "; 00585 else 00586 std::cout << 0 << " "; 00587 } 00588 std::cout << std::endl; 00589 } 00590 #endif 00591 // now get the status of neighbors; if already solved, make them 0, so not to bother 00592 // anymore on that side of tgt 00593 for( int j = 0; j < nsidesTgt; j++ ) 00594 { 00595 EntityHandle tgtNeigh = tgtNeighbors[j]; 00596 unsigned char status = 1; 00597 if( tgtNeigh == 0 ) continue; 00598 rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval ); // status 0 is unused 00599 if( 1 == status ) tgtNeighbors[j] = 0; // so will not look anymore on this side of tgt 00600 } 00601 00602 #ifdef ENABLE_DEBUG 00603 if( dbg_1 ) 00604 { 00605 std::cout << "reset sources: "; 00606 for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr ) 00607 std::cout << mb->id_from_handle( *itr ) << " "; 00608 std::cout << std::endl; 00609 } 00610 #endif 00611 EntityHandle currentSrc = srcQueue.front(); 00612 // tgt and src queues are parallel; for clarity we should have kept in the queue pairs 00613 // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with 00614 // pairs; 00615 // at every moment, the queue contains pairs of cells that intersect, and they form the 00616 // "advancing front" 00617 srcQueue.pop(); 00618 toResetSrcs.clear(); // empty the range of used srcs, will have to be set unused again, 00619 // at the end of tgt element processing 00620 toResetSrcs.insert( currentSrc ); 00621 // mb2->set_tag_data 00622 std::queue< EntityHandle > localSrc; 00623 localSrc.push( currentSrc ); 00624 #ifdef VERBOSE 00625 int countingStart = counting; 00626 #endif 00627 // will advance-front search in the neighborhood of tgt cell, until we finish processing 00628 // all 00629 // possible src cells; localSrc queue will contain all possible src cells that cover 00630 // the current tgt cell 00631 while( !localSrc.empty() ) 00632 { 00633 // 00634 EntityHandle srcT = localSrc.front(); 00635 localSrc.pop(); 00636 double P[10 * MAXEDGES], area; // 00637 int nP = 0; 00638 int nb[MAXEDGES] = { 0 }; 00639 int nr[MAXEDGES] = { 0 }; 00640 00641 int nsidesSrc; /// 00642 // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane 00643 // intersection points could include the vertices of initial elements 00644 // nb [j] = 0 means no intersection on the side j for element src (markers) 00645 // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the 00646 // tgt poly. A potential next poly in the tgt queue is the tgt poly that is 00647 // adjacent to this side 00648 rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc, 00649 nsidesTgt );MB_CHK_ERR( rval ); 00650 if( nP > 0 ) 00651 { 00652 #ifdef ENABLE_DEBUG 00653 if( dbg_1 ) 00654 { 00655 for( int k = 0; k < 3; k++ ) 00656 { 00657 std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n"; 00658 } 00659 } 00660 #endif 00661 00662 // intersection found: output P and original triangles if nP > 2 00663 EntityHandle neighbors[MAXEDGES] = { 0 }; 00664 rval = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors ); 00665 if( rval != MB_SUCCESS ) 00666 { 00667 std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT ); 00668 return MB_FAILURE; 00669 } 00670 00671 // add neighbors to the localSrc queue, if they are not marked 00672 for( int nn = 0; nn < nsidesSrc; nn++ ) 00673 { 00674 EntityHandle neighbor = neighbors[nn]; 00675 if( neighbor > 0 && nb[nn] > 0 ) // advance across src boundary nn 00676 { 00677 if( toResetSrcs.find( neighbor ) == toResetSrcs.end() ) 00678 { 00679 localSrc.push( neighbor ); 00680 #ifdef ENABLE_DEBUG 00681 if( dbg_1 ) 00682 { 00683 std::cout << " local src elem " << mb->id_from_handle( neighbor ) 00684 << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n"; 00685 mb->list_entities( &neighbor, 1 ); 00686 } 00687 #endif 00688 toResetSrcs.insert( neighbor ); 00689 } 00690 } 00691 } 00692 // n(find(nc>0))=ac; % ac is starting candidate for neighbor 00693 for( int nn = 0; nn < nsidesTgt; nn++ ) 00694 { 00695 if( nr[nn] > 0 && tgtNeighbors[nn] > 0 ) 00696 nextSrc[nn].insert( srcT ); // potential src cell that can intersect 00697 // the tgt neighbor nn 00698 } 00699 if( nP > 1 ) 00700 { // this will also construct triangles/polygons in the new mesh, if needed 00701 rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval ); 00702 } 00703 00704 recoveredArea += area; 00705 } 00706 #ifdef ENABLE_DEBUG 00707 else if( dbg_1 ) 00708 { 00709 std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " " 00710 << mb->id_from_handle( srcT ) << "\n"; 00711 } 00712 #endif 00713 } // end while (!localSrc.empty()) 00714 recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell; // replace now with recovery fraction 00715 #if defined( ENABLE_DEBUG ) || defined( VERBOSE ) 00716 if( fabs( recoveredArea ) > epsilon_1 ) 00717 { 00718 #ifdef VERBOSE 00719 std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell ) 00720 << " fraction error recovery:" << recoveredArea 00721 << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart 00722 << "\n"; 00723 #endif 00724 } 00725 #endif 00726 // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival 00727 // mesh) 00728 rs22.erase( currentTgt ); 00729 // also, look at its neighbors, and add to the seeds a next one 00730 00731 for( int j = 0; j < nsidesTgt; j++ ) 00732 { 00733 EntityHandle tgtNeigh = tgtNeighbors[j]; 00734 if( tgtNeigh == 0 || nextSrc[j].size() == 0 ) // if tgt is bigger than src, there could be no src 00735 // to advance on that side 00736 continue; 00737 int nsidesTgt2 = 0; 00738 setup_tgt_cell( tgtNeigh, 00739 nsidesTgt2 ); // find possible intersection with src cell from nextSrc 00740 for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit ) 00741 { 00742 EntityHandle nextB = *nit; 00743 // we identified tgt quad n[j] as possibly intersecting with neighbor j of the 00744 // src quad 00745 double P[10 * MAXEDGES], area; // 00746 int nP = 0; 00747 int nb[MAXEDGES] = { 0 }; 00748 int nr[MAXEDGES] = { 0 }; 00749 00750 int nsidesSrc; /// 00751 rval = computeIntersectionBetweenTgtAndSrc( 00752 /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval ); 00753 if( area > 0 ) 00754 { 00755 tgtQueue.push( tgtNeigh ); 00756 srcQueue.push( nextB ); 00757 #ifdef ENABLE_DEBUG 00758 if( dbg_1 ) 00759 std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " " 00760 << mb->id_from_handle( nextB ) << std::endl; 00761 #endif 00762 rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval ); 00763 break; // so we are done with this side of tgt, we have found a proper next 00764 // seed 00765 } 00766 } 00767 } 00768 00769 } // end while (!tgtQueue.empty()) 00770 } 00771 #ifdef ENABLE_DEBUG 00772 if( dbg_1 ) 00773 { 00774 for( int k = 0; k < 6; k++ ) 00775 mout_1[k].close(); 00776 } 00777 #endif 00778 // before cleaning up , we need to settle the position of the intersection points 00779 // on the boundary edges 00780 // this needs to be collective, so we should maybe wait something 00781 #ifdef MOAB_HAVE_MPI 00782 rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" ); 00783 #endif 00784 00785 this->clean(); 00786 return MB_SUCCESS; 00787 } 00788 00789 // clean some memory allocated 00790 void Intx2Mesh::clean() 00791 { 00792 // 00793 int indx = 0; 00794 for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ ) 00795 { 00796 delete extraNodesVec[indx]; 00797 } 00798 // extraNodesMap.clear(); 00799 extraNodesVec.clear(); 00800 // also, delete some bit tags, used to mark processed tgts and srcs 00801 mb->tag_delete( TgtFlagTag ); 00802 counting = 0; // reset counting to original value 00803 } 00804 // this method will reduce number of nodes, collapse edges that are of length 0 00805 // so a polygon like 428 431 431 will become a line 428 431 00806 // or something like 428 431 431 531 -> 428 431 531 00807 void Intx2Mesh::correct_polygon( EntityHandle* nodes, int& nP ) 00808 { 00809 int i = 0; 00810 while( i < nP ) 00811 { 00812 int nextIndex = ( i + 1 ) % nP; 00813 if( nodes[i] == nodes[nextIndex] ) 00814 { 00815 #ifdef ENABLE_DEBUG 00816 // we need to reduce nP, and collapse nodes 00817 if( dbg_1 ) 00818 { 00819 std::cout << " nodes duplicated in list: "; 00820 for( int j = 0; j < nP; j++ ) 00821 std::cout << nodes[j] << " "; 00822 std::cout << "\n"; 00823 std::cout << " node " << nodes[i] << " at index " << i << " is duplicated" 00824 << "\n"; 00825 } 00826 #endif 00827 // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0, 00828 // then next thing does nothing 00829 // (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3 00830 for( int k = i; k < nP - 1; k++ ) 00831 nodes[k] = nodes[k + 1]; 00832 nP--; // decrease the number of nodes; also, decrease i, just if we may need to check 00833 // again 00834 i--; 00835 } 00836 i++; 00837 } 00838 return; 00839 } 00840 #ifdef MOAB_HAVE_MPI 00841 ErrorCode Intx2Mesh::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts ) 00842 { 00843 localEnts.clear(); 00844 ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" ); 00845 00846 rval = mb->get_connectivity( localEnts, local_verts ); 00847 int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" ); 00848 00849 assert( parcomm != NULL ); 00850 00851 // get the position of local vertices, and decide local boxes (allBoxes...) 00852 double bmin[3] = { DBL_MAX, DBL_MAX, DBL_MAX }; 00853 double bmax[3] = { -DBL_MAX, -DBL_MAX, -DBL_MAX }; 00854 00855 std::vector< double > coords( 3 * num_local_verts ); 00856 rval = mb->get_coords( local_verts, &coords[0] );ERRORR( rval, "can't get coords of vertices " ); 00857 00858 for( int i = 0; i < num_local_verts; i++ ) 00859 { 00860 for( int k = 0; k < 3; k++ ) 00861 { 00862 double val = coords[3 * i + k]; 00863 if( val < bmin[k] ) bmin[k] = val; 00864 if( val > bmax[k] ) bmax[k] = val; 00865 } 00866 } 00867 int numprocs = parcomm->proc_config().proc_size(); 00868 allBoxes.resize( 6 * numprocs ); 00869 00870 my_rank = parcomm->proc_config().proc_rank(); 00871 for( int k = 0; k < 3; k++ ) 00872 { 00873 allBoxes[6 * my_rank + k] = bmin[k]; 00874 allBoxes[6 * my_rank + 3 + k] = bmax[k]; 00875 } 00876 00877 // now communicate to get all boxes 00878 int mpi_err; 00879 #if( MPI_VERSION >= 2 ) 00880 // use "in place" option 00881 mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE, 00882 parcomm->proc_config().proc_comm() ); 00883 #else 00884 { 00885 std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() ); 00886 mpi_err = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE, 00887 parcomm->proc_config().proc_comm() ); 00888 allBoxes = allBoxes_tmp; 00889 } 00890 #endif 00891 if( MPI_SUCCESS != mpi_err ) return MB_FAILURE; 00892 00893 if( my_rank == 0 ) 00894 { 00895 std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2 00896 << " on second mesh \n"; 00897 for( int i = 0; i < numprocs; i++ ) 00898 { 00899 std::cout << "proc: " << i << " box min: " << allBoxes[6 * i] << " " << allBoxes[6 * i + 1] << " " 00900 << allBoxes[6 * i + 2] << " \n"; 00901 std::cout << " box max: " << allBoxes[6 * i + 3] << " " << allBoxes[6 * i + 4] << " " 00902 << allBoxes[6 * i + 5] << " \n"; 00903 } 00904 } 00905 00906 return MB_SUCCESS; 00907 } 00908 ErrorCode Intx2Mesh::create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set ) 00909 { 00910 // compute the bounding box on each proc 00911 assert( parcomm != NULL ); 00912 00913 localEnts.clear(); 00914 ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" ); 00915 00916 Tag dpTag = 0; 00917 std::string tag_name( "DP" ); 00918 rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );ERRORR( rval, "can't get DP tag" ); 00919 00920 EntityHandle dum = 0; 00921 Tag corrTag; 00922 rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );ERRORR( rval, "can't get CORR tag" ); 00923 // get all local verts 00924 Range local_verts; 00925 rval = mb->get_connectivity( localEnts, local_verts ); 00926 int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" ); 00927 00928 rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );ERRORR( rval, "can't build processor boxes" ); 00929 00930 std::vector< int > gids( num_local_verts ); 00931 rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" ); 00932 00933 // now see the departure points; to what boxes should we send them? 00934 std::vector< double > dep_points( 3 * num_local_verts ); 00935 rval = mb->tag_get_data( dpTag, local_verts, (void*)&dep_points[0] );ERRORR( rval, "can't get DP tag values" ); 00936 // ranges to send to each processor; will hold vertices and elements (quads?) 00937 // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances) 00938 std::map< int, Range > Rto; 00939 int numprocs = parcomm->proc_config().proc_size(); 00940 00941 for( Range::iterator eit = localEnts.begin(); eit != localEnts.end(); ++eit ) 00942 { 00943 EntityHandle q = *eit; 00944 const EntityHandle* conn4; 00945 int num_nodes; 00946 rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" ); 00947 CartVect qbmin( DBL_MAX ); 00948 CartVect qbmax( -DBL_MAX ); 00949 for( int i = 0; i < num_nodes; i++ ) 00950 { 00951 EntityHandle v = conn4[i]; 00952 size_t index = local_verts.find( v ) - local_verts.begin(); 00953 CartVect dp( &dep_points[3 * index] ); // will use constructor 00954 for( int j = 0; j < 3; j++ ) 00955 { 00956 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j]; 00957 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j]; 00958 } 00959 } 00960 for( int p = 0; p < numprocs; p++ ) 00961 { 00962 CartVect bbmin( &allBoxes[6 * p] ); 00963 CartVect bbmax( &allBoxes[6 * p + 3] ); 00964 if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) ) { Rto[p].insert( q ); } 00965 } 00966 } 00967 00968 // now, build TLv and TLq, for each p 00969 size_t numq = 0; 00970 size_t numv = 0; 00971 for( int p = 0; p < numprocs; p++ ) 00972 { 00973 if( p == (int)my_rank ) continue; // do not "send" it, because it is already here 00974 Range& range_to_P = Rto[p]; 00975 // add the vertices to it 00976 if( range_to_P.empty() ) continue; // nothing to send to proc p 00977 Range vertsToP; 00978 rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" ); 00979 numq = numq + range_to_P.size(); 00980 numv = numv + vertsToP.size(); 00981 range_to_P.merge( vertsToP ); 00982 } 00983 TupleList TLv; 00984 TupleList TLq; 00985 TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, DP points 00986 TLv.enableWriteAccess(); 00987 00988 int sizeTuple = 2 + max_edges_1; // determined earlier, for src, first mesh 00989 TLq.initialize( 2 + max_edges_1, 0, 1, 0, 00990 numq ); // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh 00991 TLq.enableWriteAccess(); 00992 #ifdef VERBOSE 00993 std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n"; 00994 #endif 00995 for( int to_proc = 0; to_proc < numprocs; to_proc++ ) 00996 { 00997 if( to_proc == (int)my_rank ) continue; 00998 Range& range_to_P = Rto[to_proc]; 00999 Range V = range_to_P.subset_by_type( MBVERTEX ); 01000 01001 for( Range::iterator it = V.begin(); it != V.end(); ++it ) 01002 { 01003 EntityHandle v = *it; 01004 unsigned int index = local_verts.find( v ) - local_verts.begin(); 01005 int n = TLv.get_n(); 01006 TLv.vi_wr[2 * n] = to_proc; // send to processor 01007 TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range 01008 TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i] 01009 TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1]; 01010 TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2]; 01011 TLv.inc_n(); 01012 } 01013 // also, prep the quad for sending ... 01014 Range Q = range_to_P.subset_by_dimension( 2 ); 01015 for( Range::iterator it = Q.begin(); it != Q.end(); ++it ) 01016 { 01017 EntityHandle q = *it; 01018 int global_id; 01019 rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" ); 01020 int n = TLq.get_n(); 01021 TLq.vi_wr[sizeTuple * n] = to_proc; // 01022 TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ... 01023 const EntityHandle* conn4; 01024 int num_nodes; 01025 rval = mb->get_connectivity( q, conn4, 01026 num_nodes ); // could be up to MAXEDGES, but it is limited by max_edges_1 01027 ERRORR( rval, "can't get connectivity for cell" ); 01028 if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" ); 01029 for( int i = 0; i < num_nodes; i++ ) 01030 { 01031 EntityHandle v = conn4[i]; 01032 unsigned int index = local_verts.find( v ) - local_verts.begin(); 01033 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index]; 01034 } 01035 for( int k = num_nodes; k < max_edges_1; k++ ) 01036 { 01037 TLq.vi_wr[sizeTuple * n + 2 + k] = 01038 0; // fill the rest of node ids with 0; we know that the node ids start from 1! 01039 } 01040 TLq.vul_wr[n] = q; // save here the entity handle, it will be communicated back 01041 // maybe we should forget about global ID 01042 TLq.inc_n(); 01043 } 01044 } 01045 01046 // now we are done populating the tuples; route them to the appropriate processors 01047 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 ); 01048 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 ); 01049 // the elements are already in localEnts; 01050 01051 // maps from global ids to new vertex and quad handles, that are added 01052 std::map< int, EntityHandle > globalID_to_handle; 01053 /*std::map<int, EntityHandle> globalID_to_eh;*/ 01054 globalID_to_eh.clear(); // need for next iteration 01055 // now, look at every TLv, and see if we have to create a vertex there or not 01056 int n = TLv.get_n(); // the size of the points received 01057 for( int i = 0; i < n; i++ ) 01058 { 01059 int globalId = TLv.vi_rd[2 * i + 1]; 01060 if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() ) 01061 { 01062 EntityHandle new_vert; 01063 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] }; 01064 rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " ); 01065 globalID_to_handle[globalId] = new_vert; 01066 } 01067 } 01068 01069 // now, all dep points should be at their place 01070 // look in the local list of q for this proc, and create all those quads and vertices if needed 01071 // it may be an overkill, but because it does not involve communication, we do it anyway 01072 Range& local = Rto[my_rank]; 01073 Range local_q = local.subset_by_dimension( 2 ); 01074 // the local should have all the vertices in local_verts 01075 for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it ) 01076 { 01077 EntityHandle q = *it; 01078 int nnodes; 01079 const EntityHandle* conn4; 01080 rval = mb->get_connectivity( q, conn4, nnodes );ERRORR( rval, "can't get connectivity of local q " ); 01081 EntityHandle new_conn[MAXEDGES]; 01082 for( int i = 0; i < nnodes; i++ ) 01083 { 01084 EntityHandle v1 = conn4[i]; 01085 unsigned int index = local_verts.find( v1 ) - local_verts.begin(); 01086 int globalId = gids[index]; 01087 if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() ) 01088 { 01089 // we need to create that vertex, at this position dep_points 01090 double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] }; 01091 EntityHandle new_vert; 01092 rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " ); 01093 globalID_to_handle[globalId] = new_vert; 01094 } 01095 new_conn[i] = globalID_to_handle[gids[index]]; 01096 } 01097 EntityHandle new_element; 01098 // 01099 EntityType entType = MBQUAD; 01100 if( nnodes > 4 ) entType = MBPOLYGON; 01101 if( nnodes < 4 ) entType = MBTRI; 01102 01103 rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new quad " ); 01104 rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" ); 01105 int gid_el; 01106 // get the global ID of the initial quad 01107 rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " ); 01108 globalID_to_eh[gid_el] = new_element; 01109 // is this redundant or not? 01110 rval = mb->tag_set_data( corrTag, &new_element, 1, &q );ERRORR( rval, "can't set corr tag on new el" ); 01111 // set the global id on new elem 01112 rval = mb->tag_set_data( gid, &new_element, 1, &gid_el );ERRORR( rval, "can't set global id tag on new el" ); 01113 } 01114 // now look at all elements received through; we do not want to duplicate them 01115 n = TLq.get_n(); // number of elements received by this processor 01116 // form the remote cells, that will be used to send the tracer info back to the originating proc 01117 remote_cells = new TupleList(); 01118 remote_cells->initialize( 2, 0, 1, 0, n ); // will not have tracer data anymore 01119 remote_cells->enableWriteAccess(); 01120 for( int i = 0; i < n; i++ ) 01121 { 01122 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1]; 01123 int from_proc = TLq.vi_wr[sizeTuple * i]; 01124 // do we already have a quad with this global ID, represented? 01125 if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() ) 01126 { 01127 // construct the conn quad 01128 EntityHandle new_conn[MAXEDGES]; 01129 int nnodes = -1; 01130 for( int j = 0; j < max_edges_1; j++ ) 01131 { 01132 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID 01133 if( vgid == 0 ) 01134 new_conn[j] = 0; 01135 else 01136 { 01137 assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() ); 01138 new_conn[j] = globalID_to_handle[vgid]; 01139 nnodes = j + 1; // nodes are at the beginning, and are variable number 01140 } 01141 } 01142 EntityHandle new_element; 01143 // 01144 EntityType entType = MBQUAD; 01145 if( nnodes > 4 ) entType = MBPOLYGON; 01146 if( nnodes < 4 ) entType = MBTRI; 01147 rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " ); 01148 globalID_to_eh[globalIdEl] = new_element; 01149 rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" ); 01150 /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q);ERRORR(rval, "can't set corr tag on new el");*/ 01151 remote_cells->vi_wr[2 * i] = from_proc; 01152 remote_cells->vi_wr[2 * i + 1] = globalIdEl; 01153 // remote_cells->vr_wr[i] = 0.; // no contribution yet sent back 01154 remote_cells->vul_wr[i] = TLq.vul_rd[i]; // this is the corresponding tgt cell (arrival) 01155 remote_cells->inc_n(); 01156 // set the global id on new elem 01157 rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set global id tag on new el" ); 01158 } 01159 } 01160 // order the remote cells tuple list, with the global id, because we will search in it 01161 // remote_cells->print("remote_cells before sorting"); 01162 moab::TupleList::buffer sort_buffer; 01163 sort_buffer.buffer_init( n ); 01164 remote_cells->sort( 1, &sort_buffer ); 01165 sort_buffer.reset(); 01166 return MB_SUCCESS; 01167 } 01168 01169 // this algorithm assumes lagr set is already created, and some elements will be coming from 01170 // other procs, and populate the covering_set 01171 // we need to keep in a tuple list the remote cells from other procs, because we need to send back 01172 // the intersection info (like area of the intx polygon, and the current concentration) maybe total 01173 // mass in that intx 01174 ErrorCode Intx2Mesh::create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set ) 01175 { 01176 EntityHandle dum = 0; 01177 01178 Tag corrTag; 01179 ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum ); 01180 // start copy from 2nd alg 01181 // compute the bounding box on each proc 01182 assert( parcomm != NULL ); 01183 if( 1 == parcomm->proc_config().proc_size() ) 01184 { 01185 covering_set = lagr_set; // nothing to communicate, it must be serial 01186 return MB_SUCCESS; 01187 } 01188 01189 // get all local verts 01190 Range local_verts; 01191 rval = mb->get_connectivity( localEnts, local_verts ); 01192 int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" ); 01193 01194 std::vector< int > gids( num_local_verts ); 01195 rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" ); 01196 01197 Range localDepCells; 01198 rval = mb->get_entities_by_dimension( lagr_set, 2, localDepCells );ERRORR( rval, "can't get ents by dimension from lagr set" ); 01199 01200 // get all lagr verts (departure vertices) 01201 Range lagr_verts; 01202 rval = mb->get_connectivity( localDepCells, lagr_verts ); // they should be created in 01203 // the same order as the euler vertices 01204 int num_lagr_verts = (int)lagr_verts.size();ERRORR( rval, "can't get local lagr vertices" ); 01205 01206 // now see the departure points position; to what boxes should we send them? 01207 std::vector< double > dep_points( 3 * num_lagr_verts ); 01208 rval = mb->get_coords( lagr_verts, &dep_points[0] );ERRORR( rval, "can't get departure points position" ); 01209 // ranges to send to each processor; will hold vertices and elements (quads?) 01210 // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances) 01211 std::map< int, Range > Rto; 01212 int numprocs = parcomm->proc_config().proc_size(); 01213 01214 for( Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit ) 01215 { 01216 EntityHandle q = *eit; 01217 const EntityHandle* conn4; 01218 int num_nodes; 01219 rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" ); 01220 CartVect qbmin( DBL_MAX ); 01221 CartVect qbmax( -DBL_MAX ); 01222 for( int i = 0; i < num_nodes; i++ ) 01223 { 01224 EntityHandle v = conn4[i]; 01225 int index = lagr_verts.index( v ); 01226 assert( -1 != index ); 01227 CartVect dp( &dep_points[3 * index] ); // will use constructor 01228 for( int j = 0; j < 3; j++ ) 01229 { 01230 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j]; 01231 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j]; 01232 } 01233 } 01234 for( int p = 0; p < numprocs; p++ ) 01235 { 01236 CartVect bbmin( &allBoxes[6 * p] ); 01237 CartVect bbmax( &allBoxes[6 * p + 3] ); 01238 if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) ) { Rto[p].insert( q ); } 01239 } 01240 } 01241 01242 // now, build TLv and TLq, for each p 01243 size_t numq = 0; 01244 size_t numv = 0; 01245 for( int p = 0; p < numprocs; p++ ) 01246 { 01247 if( p == (int)my_rank ) continue; // do not "send" it, because it is already here 01248 Range& range_to_P = Rto[p]; 01249 // add the vertices to it 01250 if( range_to_P.empty() ) continue; // nothing to send to proc p 01251 Range vertsToP; 01252 rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" ); 01253 numq = numq + range_to_P.size(); 01254 numv = numv + vertsToP.size(); 01255 range_to_P.merge( vertsToP ); 01256 } 01257 TupleList TLv; 01258 TupleList TLq; 01259 TLv.initialize( 2, 0, 0, 3, numv ); // to proc, GLOBAL ID, DP points 01260 TLv.enableWriteAccess(); 01261 01262 int sizeTuple = 2 + max_edges_1; // max edges could be up to MAXEDGES :) for polygons 01263 TLq.initialize( 2 + max_edges_1, 0, 1, 0, 01264 numq ); // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v) 01265 // send also the corresponding tgt cell it will come to 01266 TLq.enableWriteAccess(); 01267 #ifdef VERBOSE 01268 std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n"; 01269 #endif 01270 01271 for( int to_proc = 0; to_proc < numprocs; to_proc++ ) 01272 { 01273 if( to_proc == (int)my_rank ) continue; 01274 Range& range_to_P = Rto[to_proc]; 01275 Range V = range_to_P.subset_by_type( MBVERTEX ); 01276 01277 for( Range::iterator it = V.begin(); it != V.end(); ++it ) 01278 { 01279 EntityHandle v = *it; 01280 int index = lagr_verts.index( v ); // will be the same index as the corresponding vertex in euler verts 01281 assert( -1 != index ); 01282 int n = TLv.get_n(); 01283 TLv.vi_wr[2 * n] = to_proc; // send to processor 01284 TLv.vi_wr[2 * n + 1] = gids[index]; // global id needs index in the local_verts range 01285 TLv.vr_wr[3 * n] = dep_points[3 * index]; // departure position, of the node local_verts[i] 01286 TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1]; 01287 TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2]; 01288 TLv.inc_n(); 01289 } 01290 // also, prep the 2d cells for sending ... 01291 Range Q = range_to_P.subset_by_dimension( 2 ); 01292 for( Range::iterator it = Q.begin(); it != Q.end(); ++it ) 01293 { 01294 EntityHandle q = *it; // this is a src cell 01295 int global_id; 01296 rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" ); 01297 int n = TLq.get_n(); 01298 TLq.vi_wr[sizeTuple * n] = to_proc; // 01299 TLq.vi_wr[sizeTuple * n + 1] = global_id; // global id of element, used to identify it ... 01300 const EntityHandle* conn4; 01301 int num_nodes; 01302 rval = mb->get_connectivity( 01303 q, conn4, num_nodes ); // could be up to 10;ERRORR( rval, "can't get connectivity for quad" ); 01304 if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" ); 01305 for( int i = 0; i < num_nodes; i++ ) 01306 { 01307 EntityHandle v = conn4[i]; 01308 int index = lagr_verts.index( v ); 01309 assert( -1 != index ); 01310 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index]; 01311 } 01312 for( int k = num_nodes; k < max_edges_1; k++ ) 01313 { 01314 TLq.vi_wr[sizeTuple * n + 2 + k] = 01315 0; // fill the rest of node ids with 0; we know that the node ids start from 1! 01316 } 01317 EntityHandle tgtCell; 01318 rval = mb->tag_get_data( corrTag, &q, 1, &tgtCell );ERRORR( rval, "can't get corresponding tgt cell for dep cell" ); 01319 TLq.vul_wr[n] = tgtCell; // this will be sent to remote_cells, to be able to come back 01320 TLq.inc_n(); 01321 } 01322 } 01323 // now we can route them to each processor 01324 // now we are done populating the tuples; route them to the appropriate processors 01325 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 ); 01326 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 ); 01327 // the elements are already in localEnts; 01328 01329 // maps from global ids to new vertex and quad handles, that are added 01330 std::map< int, EntityHandle > globalID_to_handle; 01331 // we already have vertices from lagr set; they are already in the processor, even before 01332 // receiving other verts from neighbors 01333 int k = 0; 01334 for( Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ ) 01335 { 01336 globalID_to_handle[gids[k]] = *vit; // a little bit of overkill 01337 // we do know that the global ids between euler and lagr verts are parallel 01338 } 01339 /*std::map<int, EntityHandle> globalID_to_eh;*/ // do we need this one? 01340 globalID_to_eh.clear(); 01341 // now, look at every TLv, and see if we have to create a vertex there or not 01342 int n = TLv.get_n(); // the size of the points received 01343 for( int i = 0; i < n; i++ ) 01344 { 01345 int globalId = TLv.vi_rd[2 * i + 1]; 01346 if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() ) 01347 { 01348 EntityHandle new_vert; 01349 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] }; 01350 rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " ); 01351 globalID_to_handle[globalId] = new_vert; 01352 } 01353 } 01354 01355 // now, all dep points should be at their place 01356 // look in the local list of 2d cells for this proc, and create all those cells if needed 01357 // it may be an overkill, but because it does not involve communication, we do it anyway 01358 Range& local = Rto[my_rank]; 01359 Range local_q = local.subset_by_dimension( 2 ); 01360 // the local should have all the vertices in lagr_verts 01361 for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it ) 01362 { 01363 EntityHandle q = *it; // these are from lagr cells, local 01364 int gid_el; 01365 rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " ); 01366 globalID_to_eh[gid_el] = q; // do we need this? maybe to just mark the ones on this processor 01367 // maybe a range of global cell ids is fine? 01368 } 01369 // now look at all elements received through; we do not want to duplicate them 01370 n = TLq.get_n(); // number of elements received by this processor 01371 // a cell should be received from one proc only; so why are we so worried about duplicated 01372 // elements? a vertex can be received from multiple sources, that is fine 01373 01374 remote_cells = new TupleList(); 01375 remote_cells->initialize( 2, 0, 1, 0, n ); // no tracers anymore in these tuples 01376 remote_cells->enableWriteAccess(); 01377 for( int i = 0; i < n; i++ ) 01378 { 01379 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1]; 01380 int from_proc = TLq.vi_rd[sizeTuple * i]; 01381 // do we already have a quad with this global ID, represented? 01382 if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() ) 01383 { 01384 // construct the conn quad 01385 EntityHandle new_conn[MAXEDGES]; 01386 int nnodes = -1; 01387 for( int j = 0; j < max_edges_1; j++ ) 01388 { 01389 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j]; // vertex global ID 01390 if( vgid == 0 ) 01391 new_conn[j] = 0; 01392 else 01393 { 01394 assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() ); 01395 new_conn[j] = globalID_to_handle[vgid]; 01396 nnodes = j + 1; // nodes are at the beginning, and are variable number 01397 } 01398 } 01399 EntityHandle new_element; 01400 // 01401 EntityType entType = MBQUAD; 01402 if( nnodes > 4 ) entType = MBPOLYGON; 01403 if( nnodes < 4 ) entType = MBTRI; 01404 rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " ); 01405 globalID_to_eh[globalIdEl] = new_element; 01406 local_q.insert( new_element ); 01407 rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set gid on new element " ); 01408 } 01409 remote_cells->vi_wr[2 * i] = from_proc; 01410 remote_cells->vi_wr[2 * i + 1] = globalIdEl; 01411 // remote_cells->vr_wr[i] = 0.; will have a different tuple for communication 01412 remote_cells->vul_wr[i] = TLq.vul_rd[i]; // this is the corresponding tgt cell (arrival) 01413 remote_cells->inc_n(); 01414 } 01415 // now, create a new set, covering_set 01416 rval = mb->create_meshset( MESHSET_SET, covering_set );ERRORR( rval, "can't create new mesh set " ); 01417 rval = mb->add_entities( covering_set, local_q );ERRORR( rval, "can't add entities to new mesh set " ); 01418 // order the remote cells tuple list, with the global id, because we will search in it 01419 // remote_cells->print("remote_cells before sorting"); 01420 moab::TupleList::buffer sort_buffer; 01421 sort_buffer.buffer_init( n ); 01422 remote_cells->sort( 1, &sort_buffer ); 01423 sort_buffer.reset(); 01424 return MB_SUCCESS; 01425 // end copy 01426 } 01427 01428 ErrorCode Intx2Mesh::resolve_intersection_sharing() 01429 { 01430 if( parcomm && parcomm->size() > 1 ) 01431 { 01432 /* 01433 moab::ParallelMergeMesh pm(parcomm, epsilon_1); 01434 ErrorCode rval = pm.merge(outSet, false, 2); // resolve only the output set, do not skip 01435 local merge, use dim 2 ERRORR(rval, "can't merge intersection "); 01436 */ 01437 // look at non-owned shared vertices, that could be part of original source set 01438 // they should be removed from intx set reference, because they might not have a 01439 // correspondent on the other task 01440 Range nonOwnedVerts; 01441 Range vertsInIntx; 01442 Range intxCells; 01443 ErrorCode rval = mb->get_entities_by_dimension( outSet, 2, intxCells );MB_CHK_ERR( rval ); 01444 rval = mb->get_connectivity( intxCells, vertsInIntx );MB_CHK_ERR( rval ); 01445 01446 rval = parcomm->filter_pstatus( vertsInIntx, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonOwnedVerts );MB_CHK_ERR( rval ); 01447 01448 // some of these vertices can be in original set 1, which was covered, transported; 01449 // but they should not be "shared" from the intx point of view, because they are not shared 01450 // with another task they might have come from coverage as a plain vertex, so losing the 01451 // sharing property ? 01452 01453 Range coverVerts; 01454 rval = mb->get_connectivity( rs1, coverVerts );MB_CHK_ERR( rval ); 01455 // find out those that are on the interface 01456 Range vertsCovInterface; 01457 rval = parcomm->filter_pstatus( coverVerts, PSTATUS_INTERFACE, PSTATUS_AND, -1, &vertsCovInterface );MB_CHK_ERR( rval ); 01458 // how many of these are in 01459 Range nodesToDuplicate = intersect( vertsCovInterface, nonOwnedVerts ); 01460 // first, get all cells connected to these vertices, from intxCells 01461 01462 Range connectedCells; 01463 rval = mb->get_adjacencies( nodesToDuplicate, 2, false, connectedCells, Interface::UNION );MB_CHK_ERR( rval ); 01464 // only those in intx set: 01465 connectedCells = intersect( connectedCells, intxCells ); 01466 // first duplicate vertices in question: 01467 std::map< EntityHandle, EntityHandle > duplicatedVerticesMap; 01468 for( Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ ) 01469 { 01470 EntityHandle vertex = *vit; 01471 double coords[3]; 01472 rval = mb->get_coords( &vertex, 1, coords );MB_CHK_ERR( rval ); 01473 EntityHandle newVertex; 01474 rval = mb->create_vertex( coords, newVertex );MB_CHK_ERR( rval ); 01475 duplicatedVerticesMap[vertex] = newVertex; 01476 } 01477 01478 // look now at connectedCells, and change their connectivities: 01479 for( Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ ) 01480 { 01481 EntityHandle intxCell = *eit; 01482 // replace connectivity 01483 std::vector< EntityHandle > connectivity; 01484 rval = mb->get_connectivity( &intxCell, 1, connectivity );MB_CHK_ERR( rval ); 01485 for( size_t i = 0; i < connectivity.size(); i++ ) 01486 { 01487 EntityHandle currentVertex = connectivity[i]; 01488 std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex ); 01489 if( mit != duplicatedVerticesMap.end() ) 01490 { 01491 connectivity[i] = mit->second; // replace connectivity directly 01492 } 01493 } 01494 int nnodes = (int)connectivity.size(); 01495 rval = mb->set_connectivity( intxCell, &connectivity[0], nnodes );MB_CHK_ERR( rval ); 01496 } 01497 } 01498 return MB_SUCCESS; 01499 } 01500 #endif /* MOAB_HAVE_MPI */ 01501 } /* namespace moab */