Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
mbIntxCheck.cpp
Go to the documentation of this file.
00001 /*
00002  * Usage: MOAB intersection check tool, for intersection on a sphere or in plane
00003  *
00004  * mpiexec -np n ./mbintx_check -s <source mesh> -t <target mesh> -i <intx mesh>  -o <source verif>
00005  * -q <target verif>
00006  *
00007  * after a run with mbtempest in parallel, that computes intersection, we can use this tool to
00008  * verify areas of intersection polygons, compared to areas of source, target elements
00009  *
00010  */
00011 
00012 #include <iostream>
00013 #include <cstdlib>
00014 #include <vector>
00015 #include <string>
00016 #include <sstream>
00017 #include <cassert>
00018 
00019 #include "moab/Core.hpp"
00020 #include "moab/IntxMesh/IntxUtils.hpp"
00021 #include "moab/ProgOptions.hpp"
00022 #include "moab/CpuTimer.hpp"
00023 #include "DebugOutput.hpp"
00024 
00025 #ifdef MOAB_HAVE_MPI
00026 // MPI includes
00027 #include "moab_mpi.h"
00028 #include "moab/ParallelComm.hpp"
00029 #include "MBParallelConventions.h"
00030 #endif
00031 
00032 using namespace moab;
00033 
00034 int main( int argc, char* argv[] )
00035 {
00036     std::stringstream sstr;
00037     // Default area_method = lHuiller; Options: Girard, GaussQuadrature (if TR is available)
00038     const IntxAreaUtils::AreaMethod areaMethod = IntxAreaUtils::lHuiller;
00039 
00040     int rank = 0, size = 1;
00041 #ifdef MOAB_HAVE_MPI
00042     MPI_Init( &argc, &argv );
00043     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00044     MPI_Comm_size( MPI_COMM_WORLD, &size );
00045 #endif
00046 
00047     std::string sourceFile, targetFile, intxFile;
00048     std::string source_verif( "outS.h5m" ), target_verif( "outt.h5m" );
00049     int sphere           = 1;
00050     int oldNamesParents  = 1;
00051     double areaErrSource = -1;
00052     double areaErrTarget = -1;
00053     ProgOptions opts;
00054 
00055     opts.addOpt< std::string >( "source,s", "source file ", &sourceFile );
00056     opts.addOpt< std::string >( "target,t", "target file ", &targetFile );
00057     opts.addOpt< std::string >( "intersection,i", "intersection file ", &intxFile );
00058     opts.addOpt< std::string >( "verif_source,v", "output source verification ", &source_verif );
00059     opts.addOpt< std::string >( "verif_target,w", "output target verification ", &target_verif );
00060     opts.addOpt< double >( "threshold_source,m", "error source threshold ", &areaErrSource );
00061     opts.addOpt< double >( "threshold_target,q", "error target threshold ", &areaErrTarget );
00062 
00063     opts.addOpt< int >( "sphere,p", "mesh on a sphere", &sphere );
00064     opts.addOpt< int >( "old_convention,n", "old names for parent tags", &oldNamesParents );
00065 
00066     opts.parseCommandLine( argc, argv );
00067     // load meshes in parallel if needed
00068     std::string opts_read = ( size == 1 ? ""
00069                                         : std::string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ) +
00070                                               std::string( ";PARALLEL_RESOLVE_SHARED_ENTS" ) );
00071 
00072     // read meshes in 3 file sets
00073     ErrorCode rval;
00074     Core moab;
00075     Interface* mb = &moab;  // global
00076     EntityHandle sset, tset, ixset;
00077 
00078     // create meshsets and load files
00079     rval = mb->create_meshset( MESHSET_SET, sset );MB_CHK_ERR( rval );
00080     rval = mb->create_meshset( MESHSET_SET, tset );MB_CHK_ERR( rval );
00081     rval = mb->create_meshset( MESHSET_SET, ixset );MB_CHK_ERR( rval );
00082     if( 0 == rank ) std::cout << "Loading source file " << sourceFile << "\n";
00083     rval = mb->load_file( sourceFile.c_str(), &sset, opts_read.c_str() );MB_CHK_SET_ERR( rval, "failed reading source file" );
00084     if( 0 == rank ) std::cout << "Loading target file " << targetFile << "\n";
00085     rval = mb->load_file( targetFile.c_str(), &tset, opts_read.c_str() );MB_CHK_SET_ERR( rval, "failed reading target file" );
00086 
00087     if( 0 == rank ) std::cout << "Loading intersection file " << intxFile << "\n";
00088     rval = mb->load_file( intxFile.c_str(), &ixset, opts_read.c_str() );MB_CHK_SET_ERR( rval, "failed reading intersection file" );
00089     double R = 1.;
00090     if( sphere )
00091     {
00092         // fix radius of both meshes, to be consistent with radius 1
00093         rval = moab::IntxUtils::ScaleToRadius( mb, sset, R );MB_CHK_ERR( rval );
00094         rval = moab::IntxUtils::ScaleToRadius( mb, tset, R );MB_CHK_ERR( rval );
00095     }
00096     Range intxCells;
00097     rval = mb->get_entities_by_dimension( ixset, 2, intxCells );MB_CHK_ERR( rval );
00098 
00099     Range sourceCells;
00100     rval = mb->get_entities_by_dimension( sset, 2, sourceCells );MB_CHK_ERR( rval );
00101 
00102     Range targetCells;
00103     rval = mb->get_entities_by_dimension( tset, 2, targetCells );MB_CHK_ERR( rval );
00104 
00105     Tag sourceParentTag;
00106     Tag targetParentTag;
00107     if( oldNamesParents )
00108     {
00109         rval = mb->tag_get_handle( "RedParent", targetParentTag );MB_CHK_SET_ERR( rval, "can't find target parent tag" );
00110         rval = mb->tag_get_handle( "BlueParent", sourceParentTag );MB_CHK_SET_ERR( rval, "can't find source parent tag" );
00111     }
00112     else
00113     {
00114         rval = mb->tag_get_handle( "TargetParent", targetParentTag );MB_CHK_SET_ERR( rval, "can't find target parent tag" );
00115         rval = mb->tag_get_handle( "SourceParent", sourceParentTag );MB_CHK_SET_ERR( rval, "can't find source parent tag" );
00116     }
00117 
00118     // error sets, for better visualization
00119     EntityHandle errorSourceSet;
00120     rval = mb->create_meshset( MESHSET_SET, errorSourceSet );MB_CHK_ERR( rval );
00121     EntityHandle errorTargetSet;
00122     rval = mb->create_meshset( MESHSET_SET, errorTargetSet );MB_CHK_ERR( rval );
00123 
00124     std::map< int, double > sourceAreas;
00125     std::map< int, double > targetAreas;
00126 
00127     std::map< int, double > sourceAreasIntx;
00128     std::map< int, double > targetAreasIntx;
00129 
00130     std::map< int, int > sourceNbIntx;
00131     std::map< int, int > targetNbIntx;
00132 
00133     std::map< int, EntityHandle > sourceMap;
00134     std::map< int, EntityHandle > targetMap;
00135 
00136     Tag gidTag = mb->globalId_tag();
00137 
00138     Tag areaTag;
00139     rval = mb->tag_get_handle( "OrigArea", 1, MB_TYPE_DOUBLE, areaTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00140 
00141     moab::IntxAreaUtils areaAdaptor( areaMethod );
00142     Range non_convex_intx_cells;
00143     for( Range::iterator eit = sourceCells.begin(); eit != sourceCells.end(); ++eit )
00144     {
00145         EntityHandle cell = *eit;
00146         const EntityHandle* verts;
00147         int num_nodes;
00148         rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_ERR( rval );
00149         if( MB_SUCCESS != rval ) return -1;
00150         std::vector< double > coords( 3 * num_nodes );
00151         // get coordinates
00152         rval = mb->get_coords( verts, num_nodes, &coords[0] );
00153         if( MB_SUCCESS != rval ) return -1;
00154         double area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, R );
00155         int sourceID;
00156         rval = mb->tag_get_data( gidTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
00157         sourceAreas[sourceID] = area;
00158         rval                  = mb->tag_set_data( areaTag, &cell, 1, &area );MB_CHK_ERR( rval );
00159         sourceMap[sourceID] = cell;
00160     }
00161     for( Range::iterator eit = targetCells.begin(); eit != targetCells.end(); ++eit )
00162     {
00163         EntityHandle cell = *eit;
00164         const EntityHandle* verts;
00165         int num_nodes;
00166         rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_ERR( rval );
00167         if( MB_SUCCESS != rval ) return -1;
00168         std::vector< double > coords( 3 * num_nodes );
00169         // get coordinates
00170         rval = mb->get_coords( verts, num_nodes, &coords[0] );
00171         if( MB_SUCCESS != rval ) return -1;
00172         double area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, R );
00173         int targetID;
00174         rval = mb->tag_get_data( gidTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
00175         targetAreas[targetID] = area;
00176         rval                  = mb->tag_set_data( areaTag, &cell, 1, &area );MB_CHK_ERR( rval );
00177         targetMap[targetID] = cell;
00178     }
00179 
00180     for( Range::iterator eit = intxCells.begin(); eit != intxCells.end(); ++eit )
00181     {
00182         EntityHandle cell = *eit;
00183         const EntityHandle* verts;
00184         int num_nodes;
00185         rval = mb->get_connectivity( cell, verts, num_nodes );MB_CHK_ERR( rval );
00186         if( MB_SUCCESS != rval ) return -1;
00187         std::vector< double > coords( 3 * num_nodes );
00188         // get coordinates
00189         rval = mb->get_coords( verts, num_nodes, &coords[0] );
00190         if( MB_SUCCESS != rval ) return -1;
00191         int check_sign   = 1;
00192         double intx_area = areaAdaptor.area_spherical_polygon( &coords[0], num_nodes, R, &check_sign );
00193 
00194         rval = mb->tag_set_data( areaTag, &cell, 1, &intx_area );
00195         ;MB_CHK_ERR( rval );
00196         int sourceID, targetID;
00197         rval = mb->tag_get_data( sourceParentTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
00198         rval = mb->tag_get_data( targetParentTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
00199 
00200         std::map< int, double >::iterator sit = sourceAreasIntx.find( sourceID );
00201         if( sit == sourceAreasIntx.end() )
00202         {
00203             sourceAreasIntx[sourceID] = intx_area;
00204             sourceNbIntx[sourceID]    = 1;
00205         }
00206         else
00207         {
00208             sourceAreasIntx[sourceID] += intx_area;
00209             sourceNbIntx[sourceID]++;
00210         }
00211 
00212         std::map< int, double >::iterator tit = targetAreasIntx.find( targetID );
00213         if( tit == targetAreasIntx.end() )
00214         {
00215             targetAreasIntx[targetID] = intx_area;
00216             targetNbIntx[targetID]    = 1;
00217         }
00218         else
00219         {
00220             targetAreasIntx[targetID] += intx_area;
00221             targetNbIntx[targetID]++;
00222         }
00223         if( intx_area < 0 )
00224         {
00225             std::cout << "negative intx area: " << intx_area << " n:" << num_nodes << " parents: " << sourceID << " "
00226                       << targetID << "\n";
00227         }
00228         if( check_sign < 0 )
00229         {
00230             non_convex_intx_cells.insert( cell );
00231             std::cout << " non-convex polygon: " << intx_area << " n:" << num_nodes << " parents: " << sourceID << " "
00232                       << targetID << "\n";
00233         }
00234     }
00235     Tag diffTag;
00236     rval = mb->tag_get_handle( "AreaDiff", 1, MB_TYPE_DOUBLE, diffTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00237 
00238     Tag countIntxCellsTag;
00239     rval = mb->tag_get_handle( "CountIntx", 1, MB_TYPE_INTEGER, countIntxCellsTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
00240 
00241     for( Range::iterator eit = sourceCells.begin(); eit != sourceCells.end(); ++eit )
00242     {
00243         EntityHandle cell = *eit;
00244 
00245         int sourceID;
00246         rval = mb->tag_get_data( gidTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
00247         double areaDiff                       = sourceAreas[sourceID];
00248         std::map< int, double >::iterator sit = sourceAreasIntx.find( sourceID );
00249         int countIntxCells                    = 0;
00250         if( sit != sourceAreasIntx.end() )
00251         {
00252             areaDiff -= sourceAreasIntx[sourceID];
00253             countIntxCells = sourceNbIntx[sourceID];
00254         }
00255         rval = mb->tag_set_data( diffTag, &cell, 1, &areaDiff );
00256         rval = mb->tag_set_data( countIntxCellsTag, &cell, 1, &countIntxCells );
00257         // add to errorSourceSet set if needed
00258         if( ( areaErrSource > 0 ) && ( fabs( areaDiff ) > areaErrSource ) )
00259         {
00260             rval = mb->add_entities( errorSourceSet, &cell, 1 );MB_CHK_ERR( rval );
00261         }
00262     }
00263     if( 0 == rank ) std::cout << "write source verification file " << source_verif << "\n";
00264     rval = mb->write_file( source_verif.c_str(), 0, 0, &sset, 1 );MB_CHK_ERR( rval );
00265     if( areaErrSource > 0 )
00266     {
00267         Range sourceErrorCells;
00268         rval = mb->get_entities_by_handle( errorSourceSet, sourceErrorCells );MB_CHK_ERR( rval );
00269         EntityHandle errorSourceIntxSet;
00270         rval = mb->create_meshset( MESHSET_SET, errorSourceIntxSet );MB_CHK_ERR( rval );
00271         if( !sourceErrorCells.empty() )
00272         {
00273             // add the intx cells that have these as source parent
00274             std::vector< int > sourceIDs;
00275             sourceIDs.resize( sourceErrorCells.size() );
00276             rval = mb->tag_get_data( gidTag, sourceErrorCells, &sourceIDs[0] );MB_CHK_SET_ERR( rval, "can't get source IDs" );
00277             std::sort( sourceIDs.begin(), sourceIDs.end() );
00278             for( Range::iterator eit = intxCells.begin(); eit != intxCells.end(); ++eit )
00279             {
00280                 EntityHandle cell = *eit;
00281                 int sourceID;
00282                 rval = mb->tag_get_data( sourceParentTag, &cell, 1, &sourceID );MB_CHK_ERR( rval );
00283                 std::vector< int >::iterator j = std::lower_bound( sourceIDs.begin(), sourceIDs.end(), sourceID );
00284                 if( ( j != sourceIDs.end() ) && ( *j == sourceID ) )
00285                 {
00286                     rval = mb->add_entities( errorSourceIntxSet, &cell, 1 );
00287                     ;MB_CHK_ERR( rval );
00288                 }
00289             }
00290             std::string filtersource     = std::string( "filt_" ) + source_verif;
00291             rval                         = mb->write_file( filtersource.c_str(), 0, 0, &errorSourceSet, 1 );
00292             std::string filterIntxsource = std::string( "filtIntx_" ) + source_verif;
00293             rval                         = mb->write_file( filterIntxsource.c_str(), 0, 0, &errorSourceIntxSet, 1 );
00294         }
00295     }
00296 
00297     for( Range::iterator eit = targetCells.begin(); eit != targetCells.end(); ++eit )
00298     {
00299         EntityHandle cell = *eit;
00300 
00301         int targetID;
00302         rval = mb->tag_get_data( gidTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
00303         double areaDiff                       = targetAreas[targetID];
00304         int countIntxCells                    = 0;
00305         std::map< int, double >::iterator sit = targetAreasIntx.find( targetID );
00306         if( sit != targetAreasIntx.end() )
00307         {
00308             areaDiff -= targetAreasIntx[targetID];
00309             countIntxCells = targetNbIntx[targetID];
00310         }
00311 
00312         rval = mb->tag_set_data( diffTag, &cell, 1, &areaDiff );
00313         rval = mb->tag_set_data( countIntxCellsTag, &cell, 1, &countIntxCells );
00314         // add to errorTargetSet set if needed
00315         if( ( areaErrTarget > 0 ) && ( fabs( areaDiff ) > areaErrTarget ) )
00316         {
00317             rval = mb->add_entities( errorTargetSet, &cell, 1 );MB_CHK_ERR( rval );
00318         }
00319     }
00320     if( 0 == rank ) std::cout << "write target verification file " << target_verif << "\n";
00321     rval = mb->write_file( target_verif.c_str(), 0, 0, &tset, 1 );MB_CHK_ERR( rval );
00322     if( areaErrTarget > 0 )
00323     {
00324         Range targetErrorCells;
00325         rval = mb->get_entities_by_handle( errorTargetSet, targetErrorCells );MB_CHK_ERR( rval );
00326         if( !targetErrorCells.empty() )
00327         {
00328             EntityHandle errorTargetIntxSet;
00329             rval = mb->create_meshset( MESHSET_SET, errorTargetIntxSet );MB_CHK_ERR( rval );
00330             // add the intx cells that have these as target parent
00331             std::vector< int > targetIDs;
00332             targetIDs.resize( targetErrorCells.size() );
00333             rval = mb->tag_get_data( gidTag, targetErrorCells, &targetIDs[0] );MB_CHK_SET_ERR( rval, "can't get target IDs" );
00334             std::sort( targetIDs.begin(), targetIDs.end() );
00335             for( Range::iterator eit = intxCells.begin(); eit != intxCells.end(); ++eit )
00336             {
00337                 EntityHandle cell = *eit;
00338                 int targetID;
00339                 rval = mb->tag_get_data( targetParentTag, &cell, 1, &targetID );MB_CHK_ERR( rval );
00340                 std::vector< int >::iterator j = std::lower_bound( targetIDs.begin(), targetIDs.end(), targetID );
00341                 if( ( j != targetIDs.end() ) && ( *j == targetID ) )
00342                 {
00343                     rval = mb->add_entities( errorTargetIntxSet, &cell, 1 );
00344                     ;MB_CHK_ERR( rval );
00345                 }
00346             }
00347             std::string filterTarget     = std::string( "filt_" ) + target_verif;
00348             rval                         = mb->write_file( filterTarget.c_str(), 0, 0, &errorTargetSet, 1 );
00349             std::string filterIntxtarget = std::string( "filtIntx_" ) + target_verif;
00350             rval                         = mb->write_file( filterIntxtarget.c_str(), 0, 0, &errorTargetIntxSet, 1 );
00351         }
00352     }
00353 
00354     if( !non_convex_intx_cells.empty() )
00355     {
00356 
00357         sourceCells.clear();
00358         targetCells.clear();
00359         for( Range::iterator it = non_convex_intx_cells.begin(); it != non_convex_intx_cells.end(); it++ )
00360         {
00361             EntityHandle cellIntx = *it;
00362             int sourceID, targetID;
00363             rval = mb->tag_get_data( sourceParentTag, &cellIntx, 1, &sourceID );MB_CHK_ERR( rval );
00364             rval = mb->tag_get_data( targetParentTag, &cellIntx, 1, &targetID );MB_CHK_ERR( rval );
00365             sourceCells.insert( sourceMap[sourceID] );
00366             targetCells.insert( targetMap[targetID] );
00367         }
00368         EntityHandle nonConvexSet;
00369         rval = mb->create_meshset( MESHSET_SET, nonConvexSet );MB_CHK_ERR( rval );
00370         rval = mb->add_entities( nonConvexSet, non_convex_intx_cells );
00371         rval = mb->write_file( "nonConvex.h5m", 0, 0, &nonConvexSet, 1 );MB_CHK_ERR( rval );
00372 
00373         EntityHandle sSet;
00374         rval = mb->create_meshset( MESHSET_SET, sSet );MB_CHK_ERR( rval );
00375         rval = mb->add_entities( sSet, sourceCells );
00376         rval = mb->write_file( "nonConvexSource.h5m", 0, 0, &sSet, 1 );MB_CHK_ERR( rval );
00377         EntityHandle tSet;
00378         rval = mb->create_meshset( MESHSET_SET, tSet );MB_CHK_ERR( rval );
00379         rval = mb->add_entities( tSet, targetCells );
00380         rval = mb->write_file( "nonConvexTarget.h5m", 0, 0, &tSet, 1 );MB_CHK_ERR( rval );
00381         rval = mb->add_entities( nonConvexSet, sourceCells );
00382         rval = mb->add_entities( nonConvexSet, targetCells );
00383         rval = mb->write_file( "nonConvexAll.h5m", 0, 0, &nonConvexSet, 1 );MB_CHK_ERR( rval );
00384     }
00385     return 0;
00386 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines