8 #define _USE_MATH_DEFINES
28 :
Intx2Mesh( mbimpl ), areaMethod( amethod ), plane( 0 ), Rsrc( 0.0 ), Rdest( 0.0 )
55 for(
int i = 1; i < nsTgt; i++ )
57 middle = 1. / nsTgt * middle;
60 for(
int j = 0; j < nsTgt; j++ )
67 for(
int j = 1; j < nsTgt - 1; j++ )
86 bool check_boxes_first )
103 if( check_boxes_first )
113 if( !overlap3d && (
plane != planeb ) )
117 if( !overlap3d &&
plane == planeb )
119 for(
int j = 0; j < nsBlue; j++ )
132 for(
int j = 0; j < nsTgt; j++ )
137 for(
int j = 0; j < nsBlue; j++ )
146 for(
int j = 0; j < nsBlue; j++ )
154 std::cout <<
"gnomonic plane: " <<
plane <<
"\n";
155 std::cout <<
" target src\n";
156 for(
int j = 0; j < nsTgt; j++ )
160 for(
int j = 0; j < nsBlue; j++ )
172 if( extraPoints >= 1 )
174 for(
int k = 0; k < nsBlue; k++ )
183 markb[( k + nsBlue - 1 ) % nsBlue] =
195 if( extraPoints >= 1 )
197 for(
int k = 0; k < nsTgt; k++ )
203 markr[( k + nsTgt - 1 ) % nsTgt] =
220 for(
int k = 1; k < nP - 1; k++ )
222 #ifdef CHECK_CONVEXITY
224 for(
int k = 0; k < nP; k++ )
226 int k1 = ( k + 1 ) % nP;
227 int k2 = ( k1 + 1 ) % nP;
229 if( orientedArea < 0 )
231 std::cout <<
" oriented area is negative: " << orientedArea <<
" k:" << k <<
" target, src:" << tgt
232 <<
" " << src <<
" \n";
255 for(
int n = 0; n < nP; n++ )
256 std::cout <<
" \t" << iP[2 * n] <<
"\t" << iP[2 * n + 1] <<
"\n";
270 std::vector< EntityHandle > foundIds;
271 foundIds.resize( nP );
272 #ifdef CHECK_CONVEXITY
277 for(
int i = 0; i < nP; i++ )
279 double* pp = &iP[2 * i];
288 for( j = 0; j < nsTgt && !found; j++ )
297 #ifdef CHECK_CONVEXITY
309 for( j = 0; j < nsBlue && !found; j++ )
319 #ifdef CHECK_CONVEXITY
334 double minArea = 1.e+38;
336 for( j = 0; j < nsTgt; j++ )
338 int j1 = ( j + 1 ) % nsTgt;
345 if( area < minArea && checkx < 2 *
epsilon_1 )
352 assert( index_min >= 0 );
362 std::cerr <<
" error in adjacent target edge: " <<
mb->
id_from_handle( adjTgtEdges[index_min] )
370 int nbExtraNodesSoFar = expts->size();
371 if( nbExtraNodesSoFar > 0 )
373 std::vector< CartVect > coords1;
374 coords1.resize( nbExtraNodesSoFar );
375 mb->
get_coords( &( *expts )[0], nbExtraNodesSoFar, &( coords1[0][0] ) );
377 for(
int k = 0; k < nbExtraNodesSoFar && !found; k++ )
380 double d2 = ( pos - coords1[k] ).
length();
384 foundIds[i] = ( *expts )[k];
385 #ifdef CHECK_CONVEXITY
399 ( *expts ).push_back( outNode );
402 foundIds[i] = outNode;
409 std::cout <<
" target quad: ";
410 for(
int j1 = 0; j1 < nsTgt; j1++ )
414 std::cout <<
" a point pp is not on a target quad " << *pp <<
" " << pp[1] <<
" target quad "
422 std::cout <<
" candidate polygon: nP" << nP <<
" plane: " <<
plane <<
"\n";
423 for(
int i1 = 0; i1 < nP; i1++ )
424 std::cout << iP[2 * i1] <<
" " << iP[2 * i1 + 1] <<
" " << foundIds[i1] <<
"\n";
431 #ifdef CHECK_CONVEXITY
462 #ifdef CHECK_CONVEXITY
464 std::vector< double > coords;
465 coords.resize( 3 * nP );
467 std::vector< CartVect > posi( nP );
470 for(
int k = 0; k < nP; k++ )
472 int k1 = ( k + 1 ) % nP;
473 int k2 = ( k1 + 1 ) % nP;
474 double orientedArea =
475 area_spherical_triangle_lHuiller( &coords[3 * k], &coords[3 * k1], &coords[3 * k2],
Rdest );
476 if( orientedArea < 0 )
478 std::cout <<
" np before 1 , 2, current " << npBefore1 <<
" " << npBefore2 <<
" " << nP <<
"\n";
479 for(
int i = 0; i < nP; i++ )
481 int nexti = ( i + 1 ) % nP;
482 double lengthEdge = ( posi[i] - posi[nexti] ).
length();
483 std::cout <<
" " << foundIds[i] <<
" edge en:" << lengthEdge <<
"\n";
485 std::cout <<
" old verts: " << oldNodes <<
" other intx:" << otherIntx <<
"\n";
487 std::cout <<
"rank:" <<
my_rank <<
" oriented area in 3d is negative: " << orientedArea <<
" k:" << k
488 <<
" target, src:" << tgt <<
" " << src <<
" \n";
496 std::cout <<
"Counting: " <<
counting <<
"\n";
497 std::cout <<
" polygon " <<
mb->
id_from_handle( polyNew ) <<
" nodes: " << nP <<
" :";
498 for(
int i1 = 0; i1 < nP; i1++ )
500 std::cout <<
" plane: " <<
plane <<
"\n";
501 std::vector< CartVect > posi( nP );
503 for(
int i1 = 0; i1 < nP; i1++ )
504 std::cout << foundIds[i1] <<
" " << posi[i1] <<
"\n";
506 std::stringstream fff;
507 fff <<
"file0" <<
counting <<
".vtk";
542 if( numTracers < 1 )
MB_CHK_SET_ERR( MB_FAILURE,
"no tracers data" );
544 std::vector< double > currentVals(
rs2.
size() * numTracers );
551 int n = remote_cells->get_n();
554 remote_cells_with_tracers =
new TupleList();
555 remote_cells_with_tracers->initialize( 2, 0, 1, numTracers,
557 remote_cells_with_tracers->enableWriteAccess();
558 for(
int i = 0; i < n; i++ )
560 remote_cells_with_tracers->vi_wr[2 * i] = remote_cells->vi_wr[2 * i];
561 remote_cells_with_tracers->vi_wr[2 * i + 1] = remote_cells->vi_wr[2 * i + 1];
563 remote_cells_with_tracers->vul_wr[i] =
564 remote_cells->vul_wr[i];
565 for(
int k = 0; k < numTracers; k++ )
566 remote_cells_with_tracers->vr_wr[numTracers * i + k] = 0;
567 remote_cells_with_tracers->inc_n();
576 std::vector< double > newValues(
rs2.
size() * numTracers,
580 double check_intx_area = 0.;
585 int srcIndex, tgtIndex;
597 check_intx_area += areap;
606 if( !remote_cells_with_tracers )
MB_CHK_SET_ERR( MB_FAILURE,
"no remote cells, failure\n" );
611 int index_in_remote = remote_cells_with_tracers->find( 1, global_id_src );
612 if( index_in_remote == -1 )
613 MB_CHK_SET_ERR( MB_FAILURE,
"can't find the global id element in remote cells\n" );
614 for(
int k = 0; k < numTracers; k++ )
615 remote_cells_with_tracers->vr_wr[index_in_remote * numTracers + k] +=
616 currentVals[numTracers * ( tgtIndex - 1 ) + k] * areap;
621 int arrTgtIndex =
rs2.
index( tgtArr );
622 if( -1 == arrTgtIndex )
MB_CHK_SET_ERR( MB_FAILURE,
"can't find the target arrival index" );
623 for(
int k = 0; k < numTracers; k++ )
624 newValues[numTracers * arrTgtIndex + k] += currentVals[( tgtIndex - 1 ) * numTracers + k] * areap;
628 MB_CHK_SET_ERR( rval,
"can't get arrival target for corresponding " );
633 if( remote_cells_with_tracers )
637 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, *remote_cells_with_tracers, 0 );
641 int n = remote_cells_with_tracers->get_n();
642 for(
int j = 0; j < n; j++ )
644 EntityHandle tgtCell = remote_cells_with_tracers->vul_rd[j];
645 int arrTgtIndex =
rs2.
index( tgtCell );
646 if( -1 == arrTgtIndex )
MB_CHK_SET_ERR( MB_FAILURE,
"can't find the target arrival index" );
647 for(
int k = 0; k < numTracers; k++ )
648 newValues[arrTgtIndex * numTracers + k] += remote_cells_with_tracers->vr_rd[j * numTracers + k];
657 std::vector< double > total_mass_local( numTracers, 0. );
658 while( iter !=
rs2.
end() )
661 double* ptrArea = (
double*)data;
662 for(
int i = 0; i < count; i++, ++iter, j++, ptrArea++ )
664 for(
int k = 0; k < numTracers; k++ )
666 total_mass_local[k] += newValues[j * numTracers + k];
667 newValues[j * numTracers + k] /= ( *ptrArea );
674 std::vector< double > total_mass( numTracers, 0. );
675 double total_intx_area = 0;
677 MPI_Reduce( &total_mass_local[0], &total_mass[0], numTracers, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
678 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
680 mpi_err = MPI_Reduce( &check_intx_area, &total_intx_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
681 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
684 for(
int k = 0; k < numTracers; k++ )
685 std::cout <<
"total mass now tracer k=" << k + 1 <<
" " << total_mass[k] <<
"\n";
686 std::cout <<
"check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI *
Rsrc *
Rsrc <<
" "
687 << total_intx_area <<
"\n";
690 if( remote_cells_with_tracers )
692 delete remote_cells_with_tracers;
693 remote_cells_with_tracers = NULL;
696 for(
int k = 0; k < numTracers; k++ )
697 std::cout <<
"total mass now tracer k=" << k + 1 <<
" " << total_mass_local[k] <<
"\n";
698 std::cout <<
"check: total intersection area: (4 * M_PI * R^2): " << 4 * M_PI *
Rsrc *
Rsrc <<
" "
699 << check_intx_area <<
"\n";
709 return Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts, gnomonic );
715 int num_local_verts = (int)local_verts.
size();
717 assert( parcomm != NULL );
719 if( num_local_verts == 0 )
723 num_local_verts = (int)local_verts.
size();
729 for(
int i = 0; i < 6; i++ )
731 gnom_box[4 * i] = gnom_box[4 * i + 1] = DBL_MAX;
732 gnom_box[4 * i + 2] = gnom_box[4 * i + 3] = -DBL_MAX;
738 std::vector< double > coords( 3 * num_local_verts );
742 std::vector< int > gnplane;
743 gnplane.resize( num_local_verts );
744 for(
int i = 0; i < num_local_verts; i++ )
746 CartVect pos( &coords[3 * i] );
771 std::vector< double > elco( 3 * nnodes );
772 std::set< int > planes;
773 for(
int i = 0; i < nnodes; i++ )
775 int ix = local_verts.
index( conn[i] );
776 planes.insert( gnplane[ix] );
777 for(
int j = 0; j < 3; j++ )
779 elco[3 * i + j] = coords[3 * ix + j];
783 for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
786 for(
int i = 0; i < nnodes; i++ )
788 CartVect pos( &elco[3 * i] );
793 for(
int k = 0; k < 2; k++ )
796 if( val < gnom_box[4 * ( pl - 1 ) + k] ) gnom_box[4 * ( pl - 1 ) + k] = val;
797 if( val > gnom_box[4 * ( pl - 1 ) + 2 + k] )
798 gnom_box[4 * ( pl - 1 ) + 2 + k] = val;
804 int numprocs = parcomm->proc_config().proc_size();
807 my_rank = parcomm->proc_config().proc_rank();
808 for(
int k = 0; k < 24; k++ )
813 #if( MPI_VERSION >= 2 )
815 mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &
allBoxes[0], 24, MPI_DOUBLE,
816 parcomm->proc_config().proc_comm() );
819 std::vector< double > allBoxes_tmp( 24 * parcomm->proc_config().proc_size() );
820 mpi_err = MPI_Allgather( &
allBoxes[24 *
my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 24, MPI_DOUBLE,
821 parcomm->proc_config().proc_comm() );
825 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
830 std::cout <<
" maximum number of vertices per cell are " <<
max_edges_1 <<
" on first mesh and " <<
max_edges_2
831 <<
" on second mesh \n";
832 for(
int i = 0; i < numprocs; i++ )
834 std::cout <<
"task: " << i <<
" \n";
835 for(
int pl = 1; pl <= 6; pl++ )
837 std::cout <<
" plane " << pl <<
" min: \t" <<
allBoxes[24 * i + 4 * ( pl - 1 )] <<
" \t"
838 <<
allBoxes[24 * i + 4 * ( pl - 1 ) + 1] <<
"\n";
839 std::cout <<
" \t max: \t" <<
allBoxes[24 * i + 4 * ( pl - 1 ) + 2] <<
" \t"
840 <<
allBoxes[24 * i + 4 * ( pl - 1 ) + 3] <<
"\n";
863 assert( parcomm != NULL );
867 if( 1 == parcomm->proc_config().proc_size() )
885 &defaultInt );
MB_CHK_SET_ERR( rval,
"can't create sending processor tag" );
891 int orig_sender = -1;
896 int size_gdofs_tag = 0;
897 std::vector< int > valsDOFs;
901 if( meshCells.size() > 0 )
903 oneCell = meshCells[0];
916 valsDOFs.resize( lenTag );
922 size_gdofs_tag = lenTag;
935 int local_int_array[2], global_int_array[2];
936 local_int_array[0] = orig_sender;
937 local_int_array[1] = size_gdofs_tag;
940 MPI_Allreduce( local_int_array, global_int_array, 2, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
941 if( MPI_SUCCESS != mpi_err )
return MB_FAILURE;
942 orig_sender = global_int_array[0];
943 size_gdofs_tag = global_int_array[1];
945 std::cout <<
"proc: " <<
my_rank <<
" size_gdofs_tag:" << size_gdofs_tag <<
"\n";
947 valsDOFs.resize( size_gdofs_tag );
951 int migrated_mesh = 0;
952 if( orig_sender != -1 ) migrated_mesh = 1;
958 size_t num_mesh_verts = mesh_verts.size();
961 std::vector< double > coords_mesh( 3 * num_mesh_verts );
965 std::vector< int > gnplane;
968 gnplane.resize( num_mesh_verts );
969 for(
size_t i = 0; i < num_mesh_verts; i++ )
971 CartVect pos( &coords_mesh[3 * i] );
978 std::vector< int > gids( num_mesh_verts );
983 std::map< int, Range > Rto;
984 int numprocs = parcomm->proc_config().proc_size();
986 for(
Range::iterator eit = meshCells.begin(); eit != meshCells.end(); ++eit )
994 std::set< int > planes;
995 std::vector< double > elco( 3 * num_nodes );
996 for(
int i = 0; i < num_nodes; i++ )
999 int index = mesh_verts.index( v );
1000 if( gnomonic ) planes.insert( gnplane[index] );
1001 for(
int j = 0; j < 3; j++ )
1003 elco[3 * i + j] = coords_mesh[3 * index + j];
1009 for( std::set< int >::iterator st = planes.begin(); st != planes.end(); st++ )
1012 double qmin[2] = { DBL_MAX, DBL_MAX };
1013 double qmax[2] = { -DBL_MAX, -DBL_MAX };
1014 for(
int i = 0; i < num_nodes; i++ )
1016 CartVect dp( &elco[3 * i] );
1021 for(
int j = 0; j < 2; j++ )
1023 if( qmin[j] > c2[j] ) qmin[j] = c2[j];
1024 if( qmax[j] < c2[j] ) qmax[j] = c2[j];
1029 for(
int p = 0; p < numprocs; p++ )
1031 double procMin1 =
allBoxes[24 * p + 4 * ( pl - 1 )];
1033 if( procMin1 >= DBL_MAX )
1035 double procMin2 =
allBoxes[24 * p + 4 * ( pl - 1 ) + 1];
1036 double procMax1 =
allBoxes[24 * p + 4 * ( pl - 1 ) + 2];
1037 double procMax2 =
allBoxes[24 * p + 4 * ( pl - 1 ) + 3];
1048 for(
int p = 0; p < numprocs; p++ )
1051 bool insert =
false;
1052 for(
int i = 0; i < num_nodes; i++ )
1054 if( box.contains_point( &elco[3 * i],
box_error ) )
1060 if( insert ) Rto[p].insert( q );
1075 for(
int p = 0; p < numprocs; p++ )
1077 if( p == (
int)
my_rank )
continue;
1078 Range& range_to_P = Rto[p];
1080 if( range_to_P.empty() )
continue;
1082 std::cout <<
" proc : " <<
my_rank <<
" to proc " << p <<
" send " << range_to_P.size() <<
" cells "
1083 <<
" psize: " << range_to_P.psize() <<
"\n";
1087 numq = numq + range_to_P.size();
1088 numv = numv + vertsToP.size();
1089 range_to_P.merge( vertsToP );
1094 TLv.initialize( 2, 0, 0, 3, numv );
1095 TLv.enableWriteAccess();
1100 TLq.initialize( sizeTuple, 0, 0, 0,
1108 TLq.enableWriteAccess();
1110 std::cout <<
"from proc " <<
my_rank <<
" send " << numv <<
" vertices and " << numq <<
" elements\n";
1113 for(
int to_proc = 0; to_proc < numprocs; to_proc++ )
1115 if( to_proc == (
int)
my_rank )
continue;
1116 Range& range_to_P = Rto[to_proc];
1117 Range V = range_to_P.subset_by_type(
MBVERTEX );
1122 int index = mesh_verts.index( v );
1123 assert( -1 != index );
1124 int n = TLv.get_n();
1125 TLv.vi_wr[2 * n] = to_proc;
1126 TLv.vi_wr[2 * n + 1] = gids[index];
1127 TLv.vr_wr[3 * n] = coords_mesh[3 * index];
1128 TLv.vr_wr[3 * n + 1] = coords_mesh[3 * index + 1];
1129 TLv.vr_wr[3 * n + 2] = coords_mesh[3 * index + 2];
1133 Range Q = range_to_P.subset_by_dimension( 2 );
1139 int n = TLq.get_n();
1140 TLq.vi_wr[sizeTuple * n] = to_proc;
1141 TLq.vi_wr[sizeTuple * n + 1] =
1150 for(
int i = 0; i < num_nodes; i++ )
1153 int index = mesh_verts.index( v );
1154 assert( -1 != index );
1155 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
1159 TLq.vi_wr[sizeTuple * n + 2 + k] =
1167 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple] = orig_sender;
1168 currentIndexIntTuple++;
1171 if( size_gdofs_tag )
1174 for(
int i = 0; i < size_gdofs_tag; i++ )
1176 TLq.vi_wr[sizeTuple * n + currentIndexIntTuple + i] =
1187 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
1188 ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
1194 std::map< int, EntityHandle > globalID_to_vertex_handle;
1200 for(
Range::iterator vit = mesh_verts.begin(); vit != mesh_verts.end(); ++vit, k++ )
1202 globalID_to_vertex_handle[gids[k]] = *vit;
1205 globalID_to_eh.clear();
1208 int n = TLv.get_n();
1209 for(
int i = 0; i < n; i++ )
1211 int globalId = TLv.vi_rd[2 * i + 1];
1212 if( globalID_to_vertex_handle.find( globalId ) ==
1213 globalID_to_vertex_handle.end() )
1217 double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
1219 globalID_to_vertex_handle[globalId] = new_vert;
1230 Range local_q = local.subset_by_dimension( 2 );
1232 for(
Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
1237 assert( gid_el >= 0 );
1238 globalID_to_eh[gid_el] = q;
1247 for(
int i = 0; i < n; i++ )
1249 int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
1260 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];
1266 assert( globalID_to_vertex_handle.find( vgid ) != globalID_to_vertex_handle.end() );
1267 new_conn[j] = globalID_to_vertex_handle[vgid];
1273 EntityType entType =
MBQUAD;
1275 if( nnodes < 4 ) entType =
MBTRI;
1278 globalID_to_eh[globalIdEl] = new_element;
1279 local_q.insert( new_element );
1284 orig_sender = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple];
1286 currentIndexIntTuple++;
1289 if( size_gdofs_tag )
1291 for(
int j = 0; j < size_gdofs_tag; j++ )
1293 valsDOFs[j] = TLq.vi_wr[sizeTuple * i + currentIndexIntTuple + j];
1295 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" );
1298 int from_proc = TLq.vi_rd[sizeTuple * i];
1305 std::cout <<
" proc " <<
my_rank <<
" add " << local_q.size() <<
" cells to covering set \n";
1312 #undef CHECK_CONVEXITY