8 #define _USE_MATH_DEFINES
24 #define CHECK_CONVEXITY
29 :
Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
56 for(
int i = 1; i < nsTgt; i++ )
58 middle = 1. / nsTgt * middle;
61 for(
int j = 0; j < nsTgt; j++ )
68 for(
int j = 1; j < nsTgt - 1; j++ )
87 bool check_boxes_first )
104 if( check_boxes_first )
114 if( !overlap3d && (
plane != planeb ) )
118 if( !overlap3d &&
plane == planeb )
120 for(
int j = 0; j < nsBlue; j++ )
133 for(
int j = 0; j < nsTgt; j++ )
138 for(
int j = 0; j < nsBlue; j++ )
147 for(
int j = 0; j < nsBlue; j++ )
155 std::cout <<
"gnomonic plane: " <<
plane <<
"\n";
156 std::cout <<
" target src\n";
157 for(
int j = 0; j < nsTgt; j++ )
161 for(
int j = 0; j < nsBlue; j++ )
173 if( extraPoints >= 1 )
175 for(
int k = 0; k < nsBlue; k++ )
184 markb[( k + nsBlue - 1 ) % nsBlue] =
196 if( extraPoints >= 1 )
198 for(
int k = 0; k < nsTgt; k++ )
204 markr[( k + nsTgt - 1 ) % nsTgt] =
221 for(
int k = 1; k < nP - 1; k++ )
223 #ifdef CHECK_CONVEXITY
225 for(
int k = 0; k < nP; k++ )
227 int k1 = ( k + 1 ) % nP;
228 int k2 = ( k1 + 1 ) % nP;
230 if( orientedArea < 0 )
232 std::cout <<
" oriented area is negative: " << orientedArea <<
" k:" << k <<
" target, src:" << tgt
233 <<
" " << src <<
" \n";
256 for(
int n = 0; n < nP; n++ )
257 std::cout <<
" \t" << iP[2 * n] <<
"\t" << iP[2 * n + 1] <<
"\n";
271 std::vector< EntityHandle > foundIds;
272 foundIds.resize( nP );
273 #ifdef CHECK_CONVEXITY
279 for(
int i = 0; i < nP; i++ )
281 double* pp = &iP[2 * i];
290 for( j = 0; j < nsTgt && !found; j++ )
299 #ifdef CHECK_CONVEXITY
311 for( j = 0; j < nsBlue && !found; j++ )
321 #ifdef CHECK_CONVEXITY
336 double minArea = 1.e+38;
338 for( j = 0; j < nsTgt; j++ )
340 int j1 = ( j + 1 ) % nsTgt;
347 if( area < minArea && checkx < 2 *
epsilon_1 )
354 assert( index_min >= 0 );
364 std::cerr <<
" error in adjacent target edge: " <<
mb->
id_from_handle( adjTgtEdges[index_min] )
372 int nbExtraNodesSoFar = expts->size();
373 if( nbExtraNodesSoFar > 0 )
375 std::vector< CartVect > coords1;
376 coords1.resize( nbExtraNodesSoFar );
377 mb->
get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
379 for(
int k = 0; k < nbExtraNodesSoFar && !found; k++ )
382 double d2 = ( pos - coords1[k] ).
length();
386 foundIds[i] = ( *expts )[k];
387 #ifdef CHECK_CONVEXITY
401 ( *expts ).push_back( outNode );
404 foundIds[i] = outNode;
411 std::cout <<
" target quad: ";
412 for(
int j1 = 0; j1 < nsTgt; j1++ )
416 std::cout <<
" a point pp is not on a target quad " << *pp <<
" " << pp[1] <<
" target quad "
424 std::cout <<
" candidate polygon: nP" << nP <<
" plane: " <<
plane <<
"\n";
425 for(
int i1 = 0; i1 < nP; i1++ )
426 std::cout << iP[2 * i1] <<
" " << iP[2 * i1 + 1] <<
" " << foundIds[i1] <<
"\n";
433 #ifdef CHECK_CONVEXITY
464 #ifdef CHECK_CONVEXITY
466 std::vector< double > coords;
467 coords.resize( 3 * nP );
469 std::vector< CartVect > posi( nP );
472 for(
int k = 0; k < nP; k++ )
474 int k1 = ( k + 1 ) % nP;
475 int k2 = ( k1 + 1 ) % nP;
476 double orientedArea =
478 if( orientedArea < 0 )
480 std::cout <<
" np before 1 , 2, current " << npBefore1 <<
" " << npBefore2 <<
" " << nP <<
"\n";
481 for(
int i = 0; i < nP; i++ )
483 int nexti = ( i + 1 ) % nP;
484 double lengthEdge = ( posi[i] - posi[nexti] ).
length();
485 std::cout <<
" " << foundIds[i] <<
" edge en:" << lengthEdge <<
"\n";
487 std::cout <<
" old verts: " << oldNodes <<
" other intx:" << otherIntx <<
"\n";
489 std::cout <<
"rank:" <<
my_rank <<
" oriented area in 3d is negative: " << orientedArea <<
" k:" << k
490 <<
" target, src:" << tgt <<
" " << src <<
" \n";
498 std::cout <<
"Counting: " <<
counting <<
"\n";
499 std::cout <<
" polygon " <<
mb->
id_from_handle( polyNew ) <<
" nodes: " << nP <<
" :";
500 for(
int i1 = 0; i1 < nP; i1++ )
502 std::cout <<
" plane: " <<
plane <<
"\n";
503 std::vector< CartVect > posi( nP );
505 for(
int i1 = 0; i1 < nP; i1++ )
506 std::cout << foundIds[i1] <<
" " << posi[i1] <<
"\n";
508 std::stringstream fff;
509 fff <<
"file0" <<
counting <<
".vtk";
544 if( numTracers < 1 )
MB_CHK_SET_ERR( MB_FAILURE,
"no tracers data" );
546 std::vector< double > currentVals(
rs2.
size() * numTracers );
553 int n = remote_cells->get_n();
556 remote_cells_with_tracers =
new TupleList();
557 remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
559 remote_cells_with_tracers->enableWriteAccess();
560 for(
int i = 0; i < n; i++ )
562 remote_cells_with_tracers->vi_wr[2 * i] = remote_cells->vi_wr[2 * i];
563 remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
565 remote_cells_with_tracers->vul_wr[i] =
566 remote_cells->vul_wr[i];
567 for(
int k = 0; k < numTracers; k++ )
568 remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0;
569 remote_cells_with_tracers->inc_n();
578 std::vector< double > newValues(
rs2.
size() * numTracers,
582 double check_intx_area = 0.;
587 int srcIndex, tgtIndex;
597 double radius =
Rsrc;
599 check_intx_area += areap;
608 if( !remote_cells_with_tracers )
MB_CHK_SET_ERR( MB_FAILURE,
"no remote cells, failure\n" );
613 int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
614 if( index_in_remote == -1 )
615 MB_CHK_SET_ERR( MB_FAILURE,
"can't find the global id element in remote cells\n" );
616 for(
int k = 0; k < numTracers; k++ )
617 remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
618 currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
623 int arrTgtIndex =
rs2.
index( tgtArr );
624 if( -1 == arrTgtIndex )
MB_CHK_SET_ERR( MB_FAILURE,
"can't find the target arrival index" );
625 for(
int k = 0; k < numTracers; k++ )
626 newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
630 MB_CHK_SET_ERR( rval,
"can't get arrival target for corresponding " );
635 if( remote_cells_with_tracers )
639 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
643 int n = remote_cells_with_tracers->get_n();
644 for(
int j = 0; j < n; j++ )
646 EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j];
647 int arrTgtIndex =
rs2.
index( tgtCell );
648 if( -1 == arrTgtIndex )
MB_CHK_SET_ERR( MB_FAILURE,
"can't find the target arrival index" );
649 for(
int k = 0; k < numTracers; k++ )
650 newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
659 std::vector< double > total_mass_local( numTracers, 0. );
660 while( iter !=
rs2.
end() )
663 double* ptrArea = (
double*)data;
664 for(
int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
666 for(
int k = 0; k < numTracers; k++ )
668 total_mass_local[k] += newValues[j * numTracers + k];
669 newValues[j * numTracers + k] /= ( *ptrArea );
676 std::vector< double > total_mass( numTracers, 0. );
677 double total_intx_area = 0;
679 MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
680 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
682 mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
683 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
686 for(
int k = 0; k < numTracers; k++ )
687 std::cout <<
"total mass now tracer k=" << k + 1 <<
" " << total_mass[k] <<
"\n";
688 std::cout <<
"check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI *
Rsrc *
Rsrc <<
" "
689 << total_intx_area <<
"\n";
692 if( remote_cells_with_tracers )
694 delete remote_cells_with_tracers;
695 remote_cells_with_tracers = NULL;
698 for(
int k = 0; k < numTracers; k++ )
699 std::cout <<
"total mass now tracer k=" << k + 1 <<
" " << total_mass_local[k] <<
"\n";
700 std::cout <<
"check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI *
Rsrc *
Rsrc <<
" "
701 << check_intx_area <<
"\n";
711 return Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts, gnomonic );
718 int num_local_verts = (int)local_verts.
size();
720 assert( parcomm != NULL );
722 if( num_local_verts == 0 )
726 num_local_verts = (int)local_verts.
size();
732 for(
int i = 0; i < 6; i++ )
734 gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX;
735 gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX;
741 std::vector< double > coords( 3 * num_local_verts );
745 std::vector< int > gnplane;
746 gnplane.resize( num_local_verts );
747 for(
int i = 0; i < num_local_verts; i++ )
749 CartVect pos( &coords[3 * i] );
774 std::vector< double > elco( 3 * nnodes );
775 std::set< int > planes;
776 for(
int i = 0; i < nnodes; i++ )
778 int ix = local_verts.
index( conn[i] );
779 planes.insert( gnplane[ix] );
780 for(
int j = 0; j < 3; j++ )
782 elco[3 * i + j] = coords[3 * ix + j];
786 for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
789 for(
int i = 0; i < nnodes; i++ )
791 CartVect pos( &elco[3 * i] );
796 for(
int k = 0; k < 2; k++ )
799 if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val;
800 if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] )
801 gnom_box[4 * ( pl - 1 ) + 2 + k] = val;
807 int numprocs = parcomm->proc_config().proc_size();
810 my_rank = parcomm->proc_config().proc_rank();
811 for(
int k = 0; k < 24; k++ )
816 #if( MPI_VERSION >= 2 )
818 mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &
allBoxes[0], 24, MPI_DOUBLE,
819 parcomm->proc_config().proc_comm() );
822 std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() );
823 mpi_err = MPI_Allgather( &
allBoxes[24 *
my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE,
824 parcomm->proc_config().proc_comm() );
828 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
833 std::cout <<
" maximum number of vertices per cell are " <<
max_edges_1 <<
" on first mesh and " <<
max_edges_2
834 <<
" on second mesh \n";
835 for(
int i = 0; i < numprocs; i++ )
837 std::cout <<
"task: " << i <<
" \n";
838 for(
int pl = 1; pl <= 6; pl++ )
840 std::cout <<
" plane " << pl <<
" min: \t" <<
allBoxes[24 * i + 4 * ( pl - 1 )] <<
" \t"
841 <<
allBoxes[24 * i + 4 * ( pl - 1 ) + 1] <<
"\n";
842 std::cout <<
" \t max: \t" <<
allBoxes[24 * i + 4 * ( pl - 1 ) + 2] <<
" \t"
843 <<
allBoxes[24 * i + 4 * ( pl - 1 ) + 3] <<
"\n";
862 int nb_ghost_layers )
870 assert( parcomm != NULL );
874 if( 1 == parcomm->proc_config().proc_size() )
892 &defaultInt );
MB_CHK_SET_ERR( rval,
"can't create sending processor tag" );
898 int orig_sender = -1;
903 int size_gdofs_tag = 0;
904 std::vector< int > valsDOFs;
908 if( meshCells.size() > 0 )
910 oneCell = meshCells[0];
923 valsDOFs.resize( lenTag );
929 size_gdofs_tag = lenTag;
942 int local_int_array[2], global_int_array[2];
943 local_int_array[0] = orig_sender;
944 local_int_array[1] = size_gdofs_tag;
947 MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
948 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
949 orig_sender = global_int_array[0];
950 size_gdofs_tag = global_int_array[1];
952 std::cout <<
"proc: " <<
my_rank <<
" size_gdofs_tag:" << size_gdofs_tag <<
"\n";
954 valsDOFs.resize( size_gdofs_tag );
958 int migrated_mesh = 0;
959 if( orig_sender != -1 ) migrated_mesh = 1;
966 size_t num_mesh_verts = mesh_verts.size();
969 std::vector< double > coords_mesh( 3 * num_mesh_verts );
973 std::vector< int > gnplane;
976 gnplane.resize( num_mesh_verts );
977 for(
size_t i = 0; i < num_mesh_verts; i++ )
979 CartVect pos( &coords_mesh[3 * i] );
986 std::vector< int > gids( num_mesh_verts );
991 std::map< int, Range > Rto;
992 int numprocs = parcomm->proc_config().proc_size();
998 if( nb_ghost_layers > 0 )
1003 double global_diag = 0;
1004 mpi_err = MPI_Allreduce( &diagonal, &global_diag, 1, MPI_DOUBLE, MPI_MAX, parcomm->proc_config().proc_comm() );
1005 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
1006 double extra_thickness = global_diag * nb_ghost_layers;
1007 if( gnomonic ) extra_thickness *= sqrt( 3. );
1010 std::cout <<
"ghost_layers:" << nb_ghost_layers <<
" max diagonal:" << global_diag
1011 <<
" extra thickness:" << extra_thickness <<
" box_error:" <<
box_error <<
"\n";
1013 for(
Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit )
1021 std::set< int > planes;
1022 std::vector< double > elco( 3 * num_nodes );
1023 for(
int i = 0; i < num_nodes; i++ )
1026 int index = mesh_verts.index( v );
1027 if( gnomonic ) planes.insert( gnplane[index] );
1028 for(
int j = 0; j < 3; j++ )
1030 elco[3 * i + j] = coords_mesh[3 * index + j];
1036 for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
1039 double qmin[2] = { DBL_MAX, DBL_MAX };
1040 double qmax[2] = { -DBL_MAX, -DBL_MAX };
1041 for(
int i = 0; i < num_nodes; i++ )
1043 CartVect dp( &elco[3 * i] );
1048 for(
int j = 0; j < 2; j++ )
1050 if( qmin[j] > c2[j] ) qmin[j] = c2[j];
1051 if( qmax[j] < c2[j] ) qmax[j] = c2[j];
1056 for(
int p = 0; p < numprocs; p++ )
1058 double procMin1 =
allBoxes[24 * p + 4 * ( pl - 1 )];
1060 if( procMin1 >= DBL_MAX )
1062 double procMin2 =
allBoxes[24 * p + 4 * ( pl - 1 ) + 1];
1063 double procMax1 =
allBoxes[24 * p + 4 * ( pl - 1 ) + 2];
1064 double procMax2 =
allBoxes[24 * p + 4 * ( pl - 1 ) + 3];
1075 for(
int p = 0; p < numprocs; p++ )
1078 bool insert =
false;
1079 for(
int i = 0; i < num_nodes; i++ )
1081 if( box.contains_point( &elco[3 * i],
box_error ) )
1087 if( insert ) Rto[p].insert( q );
1102 for(
int p = 0; p < numprocs; p++ )
1104 if( p == (
int)
my_rank )
continue;
1105 Range& range_to_P = Rto[p];
1107 if( range_to_P.empty() )
continue;
1109 std::cout <<
" proc : " <<
my_rank <<
" to proc " << p <<
" send " << range_to_P.size() <<
" cells "
1110 <<
" psize: " << range_to_P.psize() <<
"\n";
1114 numq = numq + range_to_P.size();
1115 numv = numv + vertsToP.size();
1116 range_to_P.merge( vertsToP );
1121 TLv.initialize( 2, 0, 0, 3, numv );
1122 TLv.enableWriteAccess();
1127 TLq.initialize( sizeTuple, 0, 0, 0,
1135 TLq.enableWriteAccess();
1137 std::cout <<
"from proc " <<
my_rank <<
" send " << numv <<
" vertices and " << numq <<
" elements\n";
1140 for(
int to_proc = 0; to_proc < numprocs; to_proc++ )
1142 if( to_proc == (
int)
my_rank )
continue;
1143 Range& range_to_P = Rto[to_proc];
1144 Range V = range_to_P.subset_by_type(
MBVERTEX );
1149 int index = mesh_verts.index( v );
1150 assert( -1 != index );
1151 int n = TLv.get_n();
1152 TLv.vi_wr[2 * n] = to_proc;
1153 TLv.vi_wr[2 * n + 1] = gids[index];
1154 TLv.vr_wr[3 * n] = coords_mesh[3 * index];
1155 TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1];
1156 TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2];
1160 Range Q = range_to_P.subset_by_dimension( 2 );
1166 int n = TLq.get_n();
1167 TLq.vi_wr[sizeTuple * n] = to_proc;
1168 TLq.vi_wr[sizeTuple * n + 1] =
1178 for(
int i = 0; i < num_nodes; i++ )
1181 int index = mesh_verts.index( v );
1182 assert( -1 != index );
1183 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1187 TLq.vi_wr[sizeTuple * n + 2 + k] =
1196 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = orig_sender;
1197 currentIndexIntTuple++;
1200 if( size_gdofs_tag )
1203 for(
int i = 0; i < size_gdofs_tag; i++ )
1205 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] =
1216 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1217 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1223 std::map< int, EntityHandle > globalID_to_vertex_handle;
1229 for(
Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ )
1231 globalID_to_vertex_handle[gids[k]] = *vit;
1234 globalID_to_eh.clear();
1237 int n = TLv.get_n();
1238 for(
int i = 0; i < n; i++ )
1240 int globalId = TLv.vi_rd[2 * i + 1];
1241 if( globalID_to_vertex_handle.find( globalId ) ==
1242 globalID_to_vertex_handle.end() )
1246 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1248 globalID_to_vertex_handle[globalId] = new_vert;
1259 Range local_q = local.subset_by_dimension( 2 );
1261 for(
Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1266 assert( gid_el >= 0 );
1267 globalID_to_eh[gid_el] = q;
1276 for(
int i = 0; i < n; i++ )
1278 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1283 if( globalID_to_eh.find( globalIdEl ) != globalID_to_eh.end() )
continue;
1289 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];
1295 assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() );
1296 new_conn[j] = globalID_to_vertex_handle[vgid];
1302 EntityType entType =
MBQUAD;
1304 if( nnodes < 4 ) entType =
MBTRI;
1307 globalID_to_eh[globalIdEl] = new_element;
1308 local_q.insert( new_element );
1313 orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
1315 currentIndexIntTuple++;
1318 int from_proc = TLq.vi_rd[sizeTuple * i];
1322 if( size_gdofs_tag )
1324 for(
int j = 0; j < size_gdofs_tag; j++ )
1326 valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j];
1328 rval =
mb->
tag_set_data( gdsTag, &new_element, 1, &valsDOFs[0] );
MB_CHK_SET_ERR( rval,
"can't set GLOBAL_DOFS data on coverage mesh" );
1335 std::cout <<
" proc " <<
my_rank <<
" add " << local_q.size() <<
" cells to covering set \n";
1342 #undef CHECK_CONVEXITY