25 int Intx2Mesh::dbg_1 = 0;
29 :
mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
30 countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), tgtConn( NULL ),
31 srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ), my_rank( 0 )
34 parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
37 max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
65 if( nnodes > max_edges ) max_edges = nnodes;
72 int local_max_edges = max_edges;
75 MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
76 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
96 unsigned char def_data_bit = 0;
112 std::vector< EntityHandle >* nv =
new std::vector< EntityHandle >;
119 &defaultInt );
MB_CHK_SET_ERR( rval,
"can't create positive tag" );
122 &defaultInt );
MB_CHK_SET_ERR( rval,
"can't create negative tag" );
133 std::vector< EntityHandle > zeroh(
max_edges_2, 0 );
145 while(
tgtConn[num_nodes - 2] ==
tgtConn[num_nodes - 1] && num_nodes > 3 )
149 for( i = 0; i < num_nodes; i++ )
152 tgtConn[( i + 1 ) % num_nodes] };
153 std::vector< EntityHandle > adj_entities;
155 if( rval !=
MB_SUCCESS || adj_entities.size() < 1 )
return rval;
156 zeroh[i] = adj_entities[0];
171 std::vector< EntityHandle > neighbors( max_edges );
172 std::vector< EntityHandle > zeroh( max_edges, 0 );
187 while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
190 for(
int i = 0; i < nsides; i++ )
194 v[1] = conn4[( i + 1 ) % nsides];
196 std::vector< EntityHandle > adjcells;
197 std::vector< EntityHandle > cellsInSet;
201 size_t siz = adjcells.size();
202 for(
size_t j = 0; j < siz; j++ )
203 if(
mb->
contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
204 siz = cellsInSet.size();
208 std::cout <<
"non manifold mesh, error" <<
mb->
list_entities( &( cellsInSet[0] ), cellsInSet.size() )
209 <<
"\n";
MB_CHK_SET_ERR( MB_FAILURE,
"non-manifold input mesh set" );
220 if( cell == cellsInSet[0] )
221 neighbors[i] = cellsInSet[1];
223 neighbors[i] = cellsInSet[0];
226 for(
int i = nsides; i < max_edges; i++ )
258 std::vector< EntityHandle >* nv =
new std::vector< EntityHandle >;
265 &defaultInt );
MB_CHK_SET_ERR( rval,
"can't create positive tag" );
268 &defaultInt );
MB_CHK_SET_ERR( rval,
"can't create negative tag" );
277 std::vector< EntityHandle > zeroh(
max_edges_2, 0 );
286 while(
tgtConn[num_nodes - 2] ==
tgtConn[num_nodes - 1] && num_nodes > 3 )
290 for( i = 0; i < num_nodes; i++ )
293 tgtConn[( i + 1 ) % num_nodes] };
294 std::vector< EntityHandle > adj_entities;
296 if( rval !=
MB_SUCCESS || adj_entities.size() < 1 )
return rval;
297 zeroh[i] = adj_entities[0];
312 double max_length = 0;
314 std::vector< double > coords;
321 while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
324 for(
int j = 0; j < nnodes; j++ )
326 int next = ( j + 1 ) % nnodes;
328 leng = ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
329 ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
330 ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
332 if( leng > max_length ) max_length = leng;
338 if( max_length < 1. )
341 tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
348 std::cout <<
" max edge length: " << max_length <<
" tolerance for kd tree: " <<
tolerance <<
"\n";
349 std::cout <<
" box overlap tolerance: " <<
box_error <<
"\n";
361 double recoveredArea = 0;
362 std::vector< double > positions;
363 positions.resize( nnodes * 3 );
368 for(
int k = 0; k < nnodes; k++ )
370 int ik = ( k + 1 ) % nnodes;
372 for(
int j = 0; j < 3; j++ )
374 double len2 = positions[3 * k + j] - positions[3 * ik + j];
377 av_len += sqrt( len1 );
379 if( nnodes > 0 ) av_len /= nnodes;
382 Range close_source_cells;
383 std::vector< EntityHandle > leaves;
384 for(
int i = 0; i < nnodes; i++ )
389 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
398 if( close_source_cells.
empty() )
400 std::cout <<
" there are no close source cells to target cell " << tcell <<
" id from handle "
421 recoveredArea += area;
424 recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;
430 rval = resolve_intersection_sharing();
MB_CHK_SET_ERR( rval,
"can't correct position, Intx2Mesh.cpp \n" );
445 std::stringstream ffs, fft;
446 ffs <<
"source_rank0" <<
my_rank <<
".vtk";
448 fft <<
"target_rank0" <<
my_rank <<
".vtk";
476 while( !rs22.
empty() )
478 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
481 std::cout <<
" possible not connected arrival mesh; my_rank: " <<
my_rank <<
" counting: " <<
counting
483 std::stringstream ffo;
488 bool seedFound =
false;
498 std::vector< double > positions;
499 positions.resize( nnodes * 3 );
504 Range close_source_cells;
505 std::vector< EntityHandle > leaves;
506 for(
int i = 0; i < nnodes; i++ )
512 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
544 #if defined( VERBOSE )
546 <<
" not intx with any source\n";
548 rs22.
erase( startTgt );
551 if( !seedFound )
continue;
553 std::queue< EntityHandle > srcQueue;
554 srcQueue.push( startSrc );
555 std::queue< EntityHandle > tgtQueue;
556 tgtQueue.push( startTgt );
563 unsigned char used = 1;
566 while( !tgtQueue.empty() )
577 double recoveredArea = 0;
585 std::cout <<
"Next: neighbors for current tgt ";
586 for(
int kk = 0; kk < nsidesTgt; kk++ )
588 if( tgtNeighbors[kk] > 0 )
591 std::cout << 0 <<
" ";
593 std::cout << std::endl;
598 for(
int j = 0; j < nsidesTgt; j++ )
601 unsigned char status = 1;
602 if( tgtNeigh == 0 )
continue;
604 if( 1 == status ) tgtNeighbors[j] = 0;
610 std::cout <<
"reset sources: ";
613 std::cout << std::endl;
625 toResetSrcs.
insert( currentSrc );
627 std::queue< EntityHandle > localSrc;
628 localSrc.push( currentSrc );
636 while( !localSrc.empty() )
660 for(
int k = 0; k < 3; k++ )
662 std::cout <<
" nb, nr: " << k <<
" " << nb[k] <<
" " <<
nr[k] <<
"\n";
672 std::cout <<
" can't get the neighbors for src element " <<
mb->
id_from_handle( srcT );
677 for(
int nn = 0; nn < nsidesSrc; nn++ )
680 if( neighbor > 0 && nb[nn] > 0 )
682 if( toResetSrcs.
find( neighbor ) == toResetSrcs.
end() )
684 localSrc.push( neighbor );
693 toResetSrcs.
insert( neighbor );
698 for(
int nn = 0; nn < nsidesTgt; nn++ )
700 if(
nr[nn] > 0 && tgtNeighbors[nn] > 0 )
701 nextSrc[nn].
insert( srcT );
709 recoveredArea += area;
714 std::cout <<
" tgt, src, do not intersect: " <<
mb->
id_from_handle( currentTgt ) <<
" "
719 recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;
720 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
724 std::cout <<
" tgt area: " << areaTgtCell <<
" recovered :" << recoveredArea * ( 1 + areaTgtCell )
725 <<
" fraction error recovery:" << recoveredArea
726 <<
" tgtID: " <<
mb->
id_from_handle( currentTgt ) <<
" countingStart:" << countingStart
733 rs22.
erase( currentTgt );
736 for(
int j = 0; j < nsidesTgt; j++ )
739 if( tgtNeigh == 0 || nextSrc[j].
size() == 0 )
757 tgtNeigh, nextB, P, nP, area, nb,
nr, nsidesSrc, nsidesTgt2 );
MB_CHK_ERR( rval );
760 tgtQueue.push( tgtNeigh );
761 srcQueue.push( nextB );
764 std::cout <<
"new polys pushed: src, tgt:" <<
mb->
id_from_handle( tgtNeigh ) <<
" "
779 for(
int k = 0; k < 6; k++ )
787 rval = resolve_intersection_sharing();
MB_CHK_SET_ERR( rval,
"can't correct position, Intx2Mesh.cpp \n" );
817 int nextIndex = ( i + 1 ) % nP;
818 if( nodes[i] == nodes[nextIndex] )
824 std::cout <<
" nodes duplicated in list: ";
825 for(
int j = 0; j < nP; j++ )
826 std::cout << nodes[j] <<
" ";
828 std::cout <<
" node " << nodes[i] <<
" at index " << i <<
" is duplicated"
835 for(
int k = i; k < nP - 1; k++ )
836 nodes[k] = nodes[k + 1];
850 if( gnomonic ) gnomonic =
false;
855 int num_local_verts = (int)local_verts.
size();
ERRORR( rval,
"can't get local vertices" );
857 assert( parcomm != NULL );
860 double bmin[3] = { DBL_MAX, DBL_MAX, DBL_MAX };
861 double bmax[3] = { -DBL_MAX, -DBL_MAX, -DBL_MAX };
863 std::vector< double > coords( 3 * num_local_verts );
864 rval =
mb->
get_coords( local_verts, &coords[0] );
ERRORR( rval,
"can't get coords of vertices " );
866 for(
int i = 0; i < num_local_verts; i++ )
868 for(
int k = 0; k < 3; k++ )
870 double val = coords[3 * i + k];
871 if( val < bmin[k] ) bmin[k] = val;
872 if( val > bmax[k] ) bmax[k] = val;
875 int numprocs = parcomm->proc_config().proc_size();
878 my_rank = parcomm->proc_config().proc_rank();
879 for(
int k = 0; k < 3; k++ )
887 #if( MPI_VERSION >= 2 )
889 mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &
allBoxes[0], 6, MPI_DOUBLE,
890 parcomm->proc_config().proc_comm() );
893 std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() );
894 mpi_err = MPI_Allgather( &
allBoxes[6 *
my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
895 parcomm->proc_config().proc_comm() );
899 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
904 std::cout <<
" maximum number of vertices per cell are " <<
max_edges_1 <<
" on first mesh and " <<
max_edges_2
905 <<
" on second mesh \n";
906 for(
int i = 0; i < numprocs; i++ )
908 std::cout <<
"proc: " << i <<
" box min: " <<
allBoxes[6 * i] <<
" " <<
allBoxes[6 * i + 1] <<
" "
910 std::cout <<
" box max: " <<
allBoxes[6 * i + 3] <<
" " <<
allBoxes[6 * i + 4] <<
" "
921 assert( parcomm != NULL );
927 std::string tag_name(
"DP" );
936 int num_local_verts = (int)local_verts.size();
ERRORR( rval,
"can't get local vertices" );
938 rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );
ERRORR( rval,
"can't build processor boxes" );
940 std::vector< int > gids( num_local_verts );
944 std::vector< double > dep_points( 3 * num_local_verts );
945 rval =
mb->
tag_get_data( dpTag, local_verts, (
void*)&dep_points[0] );
ERRORR( rval,
"can't get DP tag values" );
948 std::map< int, Range > Rto;
949 int numprocs = parcomm->proc_config().proc_size();
957 CartVect qbmin( DBL_MAX );
958 CartVect qbmax( -DBL_MAX );
959 for(
int i = 0; i < num_nodes; i++ )
962 size_t index = local_verts.find( v ) - local_verts.begin();
963 CartVect dp( &dep_points[3 * index] );
964 for(
int j = 0; j < 3; j++ )
966 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
967 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
970 for(
int p = 0; p < numprocs; p++ )
973 CartVect bbmax( &
allBoxes[6 * p + 3] );
984 for(
int p = 0; p < numprocs; p++ )
986 if( p == (
int)
my_rank )
continue;
987 Range& range_to_P = Rto[p];
989 if( range_to_P.empty() )
continue;
992 numq = numq + range_to_P.size();
993 numv = numv + vertsToP.size();
994 range_to_P.merge( vertsToP );
998 TLv.initialize( 2, 0, 0, 3, numv );
999 TLv.enableWriteAccess();
1004 TLq.enableWriteAccess();
1006 std::cout <<
"from proc " <<
my_rank <<
" send " << numv <<
" vertices and " << numq <<
" elements\n";
1008 for(
int to_proc = 0; to_proc < numprocs; to_proc++ )
1010 if( to_proc == (
int)
my_rank )
continue;
1011 Range& range_to_P = Rto[to_proc];
1012 Range V = range_to_P.subset_by_type(
MBVERTEX );
1017 unsigned int index = local_verts.find( v ) - local_verts.begin();
1018 int n = TLv.get_n();
1019 TLv.vi_wr[2 * n] = to_proc;
1020 TLv.vi_wr[2 * n + 1] = gids[index];
1021 TLv.vr_wr[3 * n] = dep_points[3 * index];
1022 TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
1023 TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
1027 Range Q = range_to_P.subset_by_dimension( 2 );
1033 int n = TLq.get_n();
1034 TLq.vi_wr[sizeTuple * n] = to_proc;
1035 TLq.vi_wr[sizeTuple * n + 1] = global_id;
1040 ERRORR( rval,
"can't get connectivity for cell" );
1041 if( num_nodes >
MAXEDGES )
ERRORR( MB_FAILURE,
"too many nodes in a polygon" );
1042 for(
int i = 0; i < num_nodes; i++ )
1045 unsigned int index = local_verts.find( v ) - local_verts.begin();
1046 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1050 TLq.vi_wr[sizeTuple * n + 2 + k] =
1060 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1061 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1065 std::map< int, EntityHandle > globalID_to_handle;
1067 globalID_to_eh.clear();
1069 int n = TLv.get_n();
1070 for(
int i = 0; i < n; i++ )
1072 int globalId = TLv.vi_rd[2 * i + 1];
1073 if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
1076 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1078 globalID_to_handle[globalId] = new_vert;
1086 Range local_q = local.subset_by_dimension( 2 );
1088 for(
Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1095 for(
int i = 0; i < nnodes; i++ )
1098 unsigned int index = local_verts.find( v1 ) - local_verts.begin();
1099 int globalId = gids[index];
1100 if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
1103 double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] };
1106 globalID_to_handle[globalId] = new_vert;
1108 new_conn[i] = globalID_to_handle[gids[index]];
1112 EntityType entType =
MBQUAD;
1114 if( nnodes < 4 ) entType =
MBTRI;
1116 rval =
mb->
create_element( entType, new_conn, nnodes, new_element );
ERRORR( rval,
"can't create new quad " );
1117 rval =
mb->
add_entities( covering_lagr_set, &new_element, 1 );
ERRORR( rval,
"can't add new element to dep set" );
1121 globalID_to_eh[gid_el] = new_element;
1123 rval =
mb->
tag_set_data( corrTag, &new_element, 1, &q );
ERRORR( rval,
"can't set corr tag on new el" );
1130 remote_cells =
new TupleList();
1131 remote_cells->initialize( 2, 0, 1, 0, n );
1132 remote_cells->enableWriteAccess();
1133 for(
int i = 0; i < n; i++ )
1135 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1136 int from_proc = TLq.vi_wr[sizeTuple * i];
1138 if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
1145 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];
1150 assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
1151 new_conn[j] = globalID_to_handle[vgid];
1157 EntityType entType =
MBQUAD;
1159 if( nnodes < 4 ) entType =
MBTRI;
1160 rval =
mb->
create_element( entType, new_conn, nnodes, new_element );
ERRORR( rval,
"can't create new element " );
1161 globalID_to_eh[globalIdEl] = new_element;
1162 rval =
mb->
add_entities( covering_lagr_set, &new_element, 1 );
ERRORR( rval,
"can't add new element to dep set" );
1164 remote_cells->vi_wr[2 * i] = from_proc;
1165 remote_cells->vi_wr[2 * i + 1] = globalIdEl;
1167 remote_cells->vul_wr[i] = TLq.vul_rd[i];
1168 remote_cells->inc_n();
1170 rval =
mb->
tag_set_data(
gid, &new_element, 1, &globalIdEl );
ERRORR( rval,
"can't set global id tag on new el" );
1176 sort_buffer.buffer_init( n );
1177 remote_cells->sort( 1, &sort_buffer );
1178 sort_buffer.
reset();
1195 assert( parcomm != NULL );
1196 if( 1 == parcomm->proc_config().proc_size() )
1198 covering_set = lagr_set;
1205 int num_local_verts = (int)local_verts.size();
ERRORR( rval,
"can't get local vertices" );
1207 std::vector< int > gids( num_local_verts );
1210 Range localDepCells;
1217 int num_lagr_verts = (int)lagr_verts.size();
ERRORR( rval,
"can't get local lagr vertices" );
1220 std::vector< double > dep_points( 3 * num_lagr_verts );
1221 rval =
mb->
get_coords( lagr_verts, &dep_points[0] );
ERRORR( rval,
"can't get departure points position" );
1224 std::map< int, Range > Rto;
1225 int numprocs = parcomm->proc_config().proc_size();
1227 for(
Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit )
1233 CartVect qbmin( DBL_MAX );
1234 CartVect qbmax( -DBL_MAX );
1235 for(
int i = 0; i < num_nodes; i++ )
1238 int index = lagr_verts.index( v );
1239 assert( -1 != index );
1240 CartVect dp( &dep_points[3 * index] );
1241 for(
int j = 0; j < 3; j++ )
1243 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
1244 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
1247 for(
int p = 0; p < numprocs; p++ )
1249 CartVect bbmin( &
allBoxes[6 * p] );
1250 CartVect bbmax( &
allBoxes[6 * p + 3] );
1261 for(
int p = 0; p < numprocs; p++ )
1263 if( p == (
int)
my_rank )
continue;
1264 Range& range_to_P = Rto[p];
1266 if( range_to_P.empty() )
continue;
1269 numq = numq + range_to_P.size();
1270 numv = numv + vertsToP.size();
1271 range_to_P.merge( vertsToP );
1275 TLv.initialize( 2, 0, 0, 3, numv );
1276 TLv.enableWriteAccess();
1282 TLq.enableWriteAccess();
1284 std::cout <<
"from proc " <<
my_rank <<
" send " << numv <<
" vertices and " << numq <<
" elements\n";
1287 for(
int to_proc = 0; to_proc < numprocs; to_proc++ )
1289 if( to_proc == (
int)
my_rank )
continue;
1290 Range& range_to_P = Rto[to_proc];
1291 Range V = range_to_P.subset_by_type(
MBVERTEX );
1296 int index = lagr_verts.index( v );
1297 assert( -1 != index );
1298 int n = TLv.get_n();
1299 TLv.vi_wr[2 * n] = to_proc;
1300 TLv.vi_wr[2 * n + 1] = gids[index];
1301 TLv.vr_wr[3 * n] = dep_points[3 * index];
1302 TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
1303 TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
1307 Range Q = range_to_P.subset_by_dimension( 2 );
1313 int n = TLq.get_n();
1314 TLq.vi_wr[sizeTuple * n] = to_proc;
1315 TLq.vi_wr[sizeTuple * n + 1] = global_id;
1319 q, conn4, num_nodes );
1320 if( num_nodes >
MAXEDGES )
ERRORR( MB_FAILURE,
"too many nodes in a polygon" );
1321 for(
int i = 0; i < num_nodes; i++ )
1324 int index = lagr_verts.index( v );
1325 assert( -1 != index );
1326 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1330 TLq.vi_wr[sizeTuple * n + 2 + k] =
1334 rval =
mb->
tag_get_data( corrTag, &q, 1, &tgtCell );
ERRORR( rval,
"can't get corresponding tgt cell for dep cell" );
1335 TLq.vul_wr[n] = tgtCell;
1341 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1342 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1346 std::map< int, EntityHandle > globalID_to_handle;
1350 for(
Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ )
1352 globalID_to_handle[gids[k]] = *vit;
1356 globalID_to_eh.clear();
1358 int n = TLv.get_n();
1359 for(
int i = 0; i < n; i++ )
1361 int globalId = TLv.vi_rd[2 * i + 1];
1362 if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
1365 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1367 globalID_to_handle[globalId] = new_vert;
1375 Range local_q = local.subset_by_dimension( 2 );
1377 for(
Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1382 globalID_to_eh[gid_el] = q;
1390 remote_cells =
new TupleList();
1391 remote_cells->initialize( 2, 0, 1, 0, n );
1392 remote_cells->enableWriteAccess();
1393 for(
int i = 0; i < n; i++ )
1395 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1396 int from_proc = TLq.vi_rd[sizeTuple * i];
1398 if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
1405 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];
1410 assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
1411 new_conn[j] = globalID_to_handle[vgid];
1417 EntityType entType =
MBQUAD;
1419 if( nnodes < 4 ) entType =
MBTRI;
1420 rval =
mb->
create_element( entType, new_conn, nnodes, new_element );
ERRORR( rval,
"can't create new element " );
1421 globalID_to_eh[globalIdEl] = new_element;
1422 local_q.insert( new_element );
1425 remote_cells->vi_wr[2 * i] = from_proc;
1426 remote_cells->vi_wr[2 * i + 1] = globalIdEl;
1428 remote_cells->vul_wr[i] = TLq.vul_rd[i];
1429 remote_cells->inc_n();
1433 rval =
mb->
add_entities( covering_set, local_q );
ERRORR( rval,
"can't add entities to new mesh set " );
1437 sort_buffer.buffer_init( n );
1438 remote_cells->sort( 1, &sort_buffer );
1439 sort_buffer.
reset();
1444 ErrorCode Intx2Mesh::resolve_intersection_sharing()
1446 if( parcomm && parcomm->size() > 1 )
1456 Range nonOwnedVerts;
1472 Range vertsCovInterface;
1475 Range nodesToDuplicate =
intersect( vertsCovInterface, nonOwnedVerts );
1478 Range connectedCells;
1481 connectedCells =
intersect( connectedCells, intxCells );
1483 std::map< EntityHandle, EntityHandle > duplicatedVerticesMap;
1484 for(
Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ )
1491 duplicatedVerticesMap[
vertex] = newVertex;
1495 for(
Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ )
1499 std::vector< EntityHandle > connectivity;
1501 for(
size_t i = 0; i < connectivity.size(); i++ )
1504 std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex );
1505 if( mit != duplicatedVerticesMap.end() )
1507 connectivity[i] = mit->second;
1510 int nnodes = (int)connectivity.size();