25 int main(
int argc,
char* argv[] )
28 std::string firstModel, secondModel, outputFile;
30 firstModel = TestDir +
"unittest/mbcslam/lagrangeHomme.vtk";
31 secondModel = TestDir +
"unittest/mbcslam/eulerHomme.vtk";
34 opts.
addOpt< std::string >(
"first,t",
"first mesh filename (source)", &firstModel );
35 opts.
addOpt< std::string >(
"second,m",
"second mesh filename (target)", &secondModel );
36 opts.
addOpt< std::string >(
"outputFile,o",
"output intersection file", &outputFile );
39 double epsrel = 1.e-12;
40 double boxeps = 1.e-4;
41 outputFile =
"intx.h5m";
42 opts.
addOpt<
double >(
"radius,R",
"radius for model intx", &
R );
43 opts.
addOpt<
double >(
"epsilon,e",
"relative error in intx", &epsrel );
44 opts.
addOpt<
double >(
"boxerror,b",
"relative error for box boundaries", &boxeps );
46 int output_fraction = 0;
47 int write_files_rank = 0;
50 opts.
addOpt<
int >(
"outputFraction,f",
"output fraction of areas", &output_fraction );
51 opts.
addOpt<
int >(
"writeFiles,w",
"write files of interest", &write_files_rank );
52 opts.
addOpt<
int >(
"kdtreeOption,k",
"use kd tree for intersection", &brute_force );
57 MPI_Init( &argc, &argv );
64 std::string optsRead = (
size == 1 ?
""
65 : std::string(
"PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
66 std::string(
";PARALLEL_RESOLVE_SHARED_ENTS" ) );
78 if( 0 ==
rank ) std::cout <<
"Loading mesh file " << firstModel <<
"\n";
80 if( 0 ==
rank ) std::cout <<
"Loading mesh file " << secondModel <<
"\n";
85 std::cout <<
"Radius: " <<
R <<
"\n";
86 std::cout <<
"relative eps: " << epsrel <<
"\n";
87 std::cout <<
"box eps: " << boxeps <<
"\n";
88 std::cout <<
" use kd tree for intersection: " << brute_force <<
"\n";
115 worker.set_parallel_comm( pcomm );
129 rval = worker.build_processor_euler_boxes( sf2, local_verts );
MB_CHK_ERR( rval );
130 if( write_files_rank )
132 std::stringstream outf;
133 outf <<
"second_mesh" <<
rank <<
".h5m";
139 double elapsed = MPI_Wtime();
141 bool gnomonic =
true;
142 rval = worker.construct_covering_set( sf1, covering_set, gnomonic );
MB_CHK_ERR( rval );
143 elapsed = MPI_Wtime() - elapsed;
144 if( 0 ==
rank ) std::cout <<
"\nTime to communicate the mesh = " << elapsed << std::endl;
148 if( output_fraction )
158 assert( area_cov_set > 0 );
161 int num_cov_cells, num_comm_cells;
164 double fraction_area = comm_area / area_cov_set;
165 double fraction_num_cells =
166 (double)num_comm_cells / num_cov_cells;
168 double max_fraction_area, max_fraction_num_cells, min_fraction_area, min_fraction_num_cells;
169 double average_fraction_area, average_fraction_num_cells;
170 MPI_Reduce( &fraction_area, &max_fraction_area, 1, MPI_DOUBLE, MPI_MAX, 0,
MPI_COMM_WORLD );
171 MPI_Reduce( &fraction_num_cells, &max_fraction_num_cells, 1, MPI_DOUBLE, MPI_MAX, 0,
MPI_COMM_WORLD );
172 MPI_Reduce( &fraction_area, &min_fraction_area, 1, MPI_DOUBLE, MPI_MIN, 0,
MPI_COMM_WORLD );
173 MPI_Reduce( &fraction_num_cells, &min_fraction_num_cells, 1, MPI_DOUBLE, MPI_MIN, 0,
MPI_COMM_WORLD );
174 MPI_Reduce( &fraction_area, &average_fraction_area, 1, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD );
175 MPI_Reduce( &fraction_num_cells, &average_fraction_num_cells, 1, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD );
176 average_fraction_area /=
size;
177 average_fraction_num_cells /=
size;
180 std::cout <<
" fraction area: min: " << min_fraction_area <<
" max: " << max_fraction_area
181 <<
" average :" << average_fraction_area <<
" \n";
182 std::cout <<
" fraction num cells: min: " << min_fraction_num_cells
183 <<
" max: " << max_fraction_num_cells <<
" average: " << average_fraction_num_cells <<
" \n";
186 if( write_files_rank )
188 std::stringstream cof;
189 cof <<
"covering_mesh" <<
rank <<
".h5m";
197 if( 0 ==
rank ) std::cout <<
"Computing intersections ..\n";
199 double elapsed = MPI_Wtime();
210 elapsed = MPI_Wtime() - elapsed;
211 if( 0 ==
rank ) std::cout <<
"\nTime to compute the intersection between meshes = " << elapsed << std::endl;
227 if( write_files_rank )
229 std::stringstream outf;
230 outf <<
"intersect" <<
rank <<
".h5m";
231 rval =
mb->
write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
235 std::cout <<
"On rank : " <<
rank <<
" arrival area: " << arrival_area <<
" intersection area:" << intx_area
236 <<
" rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) <<
"\n";
239 #ifdef MOAB_HAVE_HDF5_PARALLEL
240 rval =
mb->
write_file( outputFile.c_str(), 0,
"PARALLEL=WRITE_PART", &outputSet, 1 );
MB_CHK_SET_ERR( rval,
"failed to write intx file" );