MOAB: Mesh Oriented datABase  (version 5.4.1)
intx_on_sphere_test.cpp
Go to the documentation of this file.
00001 /*
00002  * intx_on_sphere_test.cpp
00003  *
00004  *  Created on: Oct 3, 2012
00005  *      Author: iulian
00006  */
00007 #include <iostream>
00008 #include <sstream>
00009 #include <ctime>
00010 #include <cstdlib>
00011 #include <cstdio>
00012 #include <cstring>
00013 #include "moab/Core.hpp"
00014 #ifdef MOAB_HAVE_MPI
00015 #include "moab/ParallelComm.hpp"
00016 #endif
00017 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00018 #include "moab/IntxMesh/IntxUtils.hpp"
00019 #include "TestUtil.hpp"
00020 #include "moab/ProgOptions.hpp"
00021 #include <cmath>
00022 
00023 using namespace moab;
00024 
00025 int main( int argc, char* argv[] )
00026 {
00027 
00028     std::string firstModel, secondModel, outputFile;
00029 
00030     firstModel  = TestDir + "unittest/mbcslam/lagrangeHomme.vtk";
00031     secondModel = TestDir + "unittest/mbcslam/eulerHomme.vtk";
00032 
00033     ProgOptions opts;
00034     opts.addOpt< std::string >( "first,t", "first mesh filename (source)", &firstModel );
00035     opts.addOpt< std::string >( "second,m", "second mesh filename (target)", &secondModel );
00036     opts.addOpt< std::string >( "outputFile,o", "output intersection file", &outputFile );
00037 
00038     double R      = 1.;  // input
00039     double epsrel = 1.e-12;
00040     double boxeps = 1.e-4;
00041     outputFile    = "intx.h5m";
00042     opts.addOpt< double >( "radius,R", "radius for model intx", &R );
00043     opts.addOpt< double >( "epsilon,e", "relative error in intx", &epsrel );
00044     opts.addOpt< double >( "boxerror,b", "relative error for box boundaries", &boxeps );
00045 
00046     int output_fraction  = 0;
00047     int write_files_rank = 0;
00048     int brute_force      = 0;
00049 
00050     opts.addOpt< int >( "outputFraction,f", "output fraction of areas", &output_fraction );
00051     opts.addOpt< int >( "writeFiles,w", "write files of interest", &write_files_rank );
00052     opts.addOpt< int >( "kdtreeOption,k", "use kd tree for intersection", &brute_force );
00053 
00054     opts.parseCommandLine( argc, argv );
00055     int rank = 0, size = 1;
00056 #ifdef MOAB_HAVE_MPI
00057     MPI_Init( &argc, &argv );
00058     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00059     MPI_Comm_size( MPI_COMM_WORLD, &size );
00060 #endif
00061 
00062     // check command line arg second grid is red, arrival, first mesh is blue, departure
00063     // will will keep the
00064     std::string optsRead = ( size == 1 ? ""
00065                                        : std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00066                                              std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) );
00067 
00068     // read meshes in 2 file sets
00069     ErrorCode rval;
00070     Core moab;
00071     Interface* mb = &moab;  // global
00072     EntityHandle sf1, sf2, outputSet;
00073 
00074     // create meshsets and load files
00075 
00076     rval = mb->create_meshset( MESHSET_SET, sf1 );MB_CHK_ERR( rval );
00077     rval = mb->create_meshset( MESHSET_SET, sf2 );MB_CHK_ERR( rval );
00078     if( 0 == rank ) std::cout << "Loading mesh file " << firstModel << "\n";
00079     rval = mb->load_file( firstModel.c_str(), &sf1, optsRead.c_str() );MB_CHK_ERR( rval );
00080     if( 0 == rank ) std::cout << "Loading mesh file " << secondModel << "\n";
00081     rval = mb->load_file( secondModel.c_str(), &sf2, optsRead.c_str() );MB_CHK_ERR( rval );
00082 
00083     if( 0 == rank )
00084     {
00085         std::cout << "Radius:  " << R << "\n";
00086         std::cout << "relative eps:  " << epsrel << "\n";
00087         std::cout << "box eps:  " << boxeps << "\n";
00088         std::cout << " use kd tree for intersection: " << brute_force << "\n";
00089     }
00090     rval = mb->create_meshset( MESHSET_SET, outputSet );MB_CHK_ERR( rval );
00091 
00092     // fix radius of both meshes, to be consistent with input R
00093     rval = moab::IntxUtils::ScaleToRadius( mb, sf1, R );MB_CHK_ERR( rval );
00094     rval = moab::IntxUtils::ScaleToRadius( mb, sf2, R );MB_CHK_ERR( rval );
00095 
00096 #if 0
00097   // std::cout << "Fix orientation etc ..\n";
00098   //IntxUtils; those calls do nothing for a good mesh
00099   rval = fix_degenerate_quads(mb, sf1);MB_CHK_ERR(rval);
00100   rval = fix_degenerate_quads(mb, sf2);MB_CHK_ERR(rval);
00101 
00102   rval = positive_orientation(mb, sf1, R);MB_CHK_ERR(rval);
00103   rval = positive_orientation(mb, sf2, R);MB_CHK_ERR(rval);
00104 #endif
00105 
00106 #ifdef MOAB_HAVE_MPI
00107     ParallelComm* pcomm = ParallelComm::get_pcomm( mb, 0 );
00108 #endif
00109     Intx2MeshOnSphere worker( mb );
00110     IntxAreaUtils areaAdaptor;
00111 
00112     worker.set_error_tolerance( R * epsrel );
00113     worker.set_box_error( boxeps );
00114 #ifdef MOAB_HAVE_MPI
00115     worker.set_parallel_comm( pcomm );
00116 #endif
00117     // worker.SetEntityType(moab::MBQUAD);
00118     worker.set_radius_source_mesh( R );
00119     worker.set_radius_destination_mesh( R );
00120     // worker.enable_debug();
00121 
00122     rval = worker.FindMaxEdges( sf1, sf2 );MB_CHK_ERR( rval );
00123 
00124     EntityHandle covering_set;
00125 #ifdef MOAB_HAVE_MPI
00126     if( size > 1 )
00127     {
00128         Range local_verts;
00129         rval = worker.build_processor_euler_boxes( sf2, local_verts );MB_CHK_ERR( rval );  // output also the local_verts
00130         if( write_files_rank )
00131         {
00132             std::stringstream outf;
00133             outf << "second_mesh" << rank << ".h5m";
00134             rval = mb->write_file( outf.str().c_str(), 0, 0, &sf2, 1 );MB_CHK_ERR( rval );
00135         }
00136     }
00137     if( size > 1 )
00138     {
00139         double elapsed = MPI_Wtime();
00140         rval           = mb->create_meshset( moab::MESHSET_SET, covering_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
00141         rval = worker.construct_covering_set( sf1, covering_set );MB_CHK_ERR( rval );  // lots of communication if mesh is distributed very differently
00142         elapsed = MPI_Wtime() - elapsed;
00143         if( 0 == rank ) std::cout << "\nTime to communicate the mesh = " << elapsed << std::endl;
00144         // area fraction of the covering set that needed to be communicated from other processors
00145         // number of elements in the covering set communicated, compared to total number of elements
00146         // in the covering set
00147         if( output_fraction )
00148         {
00149             EntityHandle comm_set;  // set with elements communicated from other tasks
00150             rval = mb->create_meshset( MESHSET_SET, comm_set );MB_CHK_ERR( rval );
00151             // see how much more different is compared to sf1
00152             rval = mb->unite_meshset( comm_set, covering_set );MB_CHK_ERR( rval );  // will have to subtract from covering set, initial set
00153             // subtract
00154             rval = mb->subtract_meshset( comm_set, sf1 );MB_CHK_ERR( rval );
00155             // compute fractions
00156             double area_cov_set = areaAdaptor.area_on_sphere( mb, covering_set, R );
00157             assert( area_cov_set > 0 );
00158             double comm_area = areaAdaptor.area_on_sphere( mb, comm_set, R );
00159             // more important is actually the number of elements communicated
00160             int num_cov_cells, num_comm_cells;
00161             rval = mb->get_number_entities_by_dimension( covering_set, 2, num_cov_cells );MB_CHK_ERR( rval );
00162             rval = mb->get_number_entities_by_dimension( comm_set, 2, num_comm_cells );MB_CHK_ERR( rval );
00163             double fraction_area = comm_area / area_cov_set;
00164             double fraction_num_cells =
00165                 (double)num_comm_cells / num_cov_cells;  // determine min, max, average of these fractions
00166 
00167             double max_fraction_area, max_fraction_num_cells, min_fraction_area, min_fraction_num_cells;
00168             double average_fraction_area, average_fraction_num_cells;
00169             MPI_Reduce( &fraction_area, &max_fraction_area, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
00170             MPI_Reduce( &fraction_num_cells, &max_fraction_num_cells, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD );
00171             MPI_Reduce( &fraction_area, &min_fraction_area, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD );
00172             MPI_Reduce( &fraction_num_cells, &min_fraction_num_cells, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD );
00173             MPI_Reduce( &fraction_area, &average_fraction_area, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00174             MPI_Reduce( &fraction_num_cells, &average_fraction_num_cells, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD );
00175             average_fraction_area /= size;
00176             average_fraction_num_cells /= size;
00177             if( rank == 0 )
00178             {
00179                 std::cout << " fraction area:      min: " << min_fraction_area << " max: " << max_fraction_area
00180                           << " average :" << average_fraction_area << " \n";
00181                 std::cout << " fraction num cells: min: " << min_fraction_num_cells
00182                           << " max: " << max_fraction_num_cells << " average: " << average_fraction_num_cells << " \n";
00183             }
00184         }
00185         if( write_files_rank )
00186         {
00187             std::stringstream cof;
00188             cof << "covering_mesh" << rank << ".h5m";
00189             rval = mb->write_file( cof.str().c_str(), 0, 0, &covering_set, 1 );MB_CHK_ERR( rval );
00190         }
00191     }
00192     else
00193 #endif
00194         covering_set = sf1;
00195 
00196     if( 0 == rank ) std::cout << "Computing intersections ..\n";
00197 #ifdef MOAB_HAVE_MPI
00198     double elapsed = MPI_Wtime();
00199 #endif
00200     if( brute_force )
00201     {
00202         rval = worker.intersect_meshes_kdtree( covering_set, sf2, outputSet );MB_CHK_SET_ERR( rval, "failed to intersect meshes with slow method" );
00203     }
00204     else
00205     {
00206         rval = worker.intersect_meshes( covering_set, sf2, outputSet );MB_CHK_SET_ERR( rval, "failed to intersect meshes" );
00207     }
00208 #ifdef MOAB_HAVE_MPI
00209     elapsed = MPI_Wtime() - elapsed;
00210     if( 0 == rank ) std::cout << "\nTime to compute the intersection between meshes = " << elapsed << std::endl;
00211 #endif
00212     // the output set does not have the intx vertices on the boundary shared, so they will be
00213     // duplicated right now we write this file just for checking it looks OK
00214 
00215     // compute total area with 2 methods
00216     // double initial_area = areaAdaptor.area_on_sphere_lHuiller(mb, sf1, R);
00217     // double area_method1 = areaAdaptor.area_on_sphere_lHuiller(mb, outputSet, R);
00218     // double area_method2 = areaAdaptor.area_on_sphere(mb, outputSet, R);
00219 
00220     // std::cout << "initial area: " << initial_area << "\n";
00221     // std::cout<< " area with l'Huiller: " << area_method1 << " with Girard: " << area_method2<<
00222     // "\n"; std::cout << " relative difference areas " <<
00223     // fabs(area_method1-area_method2)/area_method1 << "\n"; std::cout << " relative error " <<
00224     // fabs(area_method1-initial_area)/area_method1 << "\n";
00225 
00226     if( write_files_rank )
00227     {
00228         std::stringstream outf;
00229         outf << "intersect" << rank << ".h5m";
00230         rval = mb->write_file( outf.str().c_str(), 0, 0, &outputSet, 1 );
00231     }
00232     double intx_area    = areaAdaptor.area_on_sphere( mb, outputSet, R );
00233     double arrival_area = areaAdaptor.area_on_sphere( mb, sf2, R );
00234     std::cout << "On rank : " << rank << " arrival area: " << arrival_area << "  intersection area:" << intx_area
00235               << " rel error: " << fabs( ( intx_area - arrival_area ) / arrival_area ) << "\n";
00236 
00237 #ifdef MOAB_HAVE_MPI
00238 #ifdef MOAB_HAVE_HDF5_PARALLEL
00239     rval = mb->write_file( outputFile.c_str(), 0, "PARALLEL=WRITE_PART", &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
00240 #else
00241     // write intx set on rank 0, in serial; we cannot write in parallel
00242     if( 0 == rank )
00243     {
00244         rval = mb->write_file( outputFile.c_str(), 0, 0, &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
00245     }
00246 #endif
00247     MPI_Finalize();
00248 #else
00249     rval = mb->write_file( outputFile.c_str(), 0, 0, &outputSet, 1 );MB_CHK_SET_ERR( rval, "failed to write intx file" );
00250 #endif
00251     return 0;
00252 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines