MOAB: Mesh Oriented datABase  (version 5.3.0)
ExtractClose.cpp
Go to the documentation of this file.
00001 /** \brief This test shows how to extract mesh from a model, based on distance.
00002  *
00003  * MOAB's It is needed to extract from large mesh files cells close to some point, where we suspect
00004  * errors It would be useful for 9Gb input file that we cannot visualize, but we have overlapped
00005  * elements.
00006  */
00007 
00008 #include <iostream>
00009 #include <cstdlib>
00010 #include <cstdio>
00011 
00012 #include "moab/Core.hpp"
00013 #include "moab/Interface.hpp"
00014 #include "moab/Range.hpp"
00015 #include "moab/ProgOptions.hpp"
00016 #include "moab/CartVect.hpp"
00017 #include "moab/IntxMesh/IntxUtils.hpp"
00018 
00019 #ifdef MOAB_HAVE_MPI
00020 #include "moab_mpi.h"
00021 #include "moab/ParallelComm.hpp"
00022 #endif
00023 
00024 using namespace moab;
00025 using namespace std;
00026 
00027 #ifdef MOAB_HAVE_HDF5
00028 string test_file_name = string( MESH_DIR ) + string( "/64bricks_512hex_256part.h5m" );
00029 #endif
00030 int main( int argc, char** argv )
00031 {
00032     ProgOptions opts;
00033 
00034     // vertex 25 in filt out  (-0.784768, 0.594952, -0.173698)
00035     double x = -0.784768, y = 0.594952, z = -0.173698;
00036     double lat = 28.6551, lon = 61.0204;  // from Charlie Zender's page
00037     // https://acme-climate.atlassian.net/wiki/spaces/ED/pages/1347092547/Assessing+and+Addressing+Regridding+Weights+in+Map-Files
00038 
00039     string inputFile = test_file_name;
00040     opts.addOpt< string >( "inFile,i", "Specify the input file name string ", &inputFile );
00041 
00042     string outFile = "out.h5m";
00043     opts.addOpt< string >( "outFile,o", "Specify the output file name string ", &outFile );
00044 
00045     opts.addOpt< double >( string( "xpos,x" ), string( "x position" ), &x );
00046     opts.addOpt< double >( string( "ypos,y" ), string( "y position" ), &y );
00047     opts.addOpt< double >( string( "zpos,z" ), string( "z position" ), &z );
00048 
00049     opts.addOpt< double >( string( "latitude,t" ), string( "north latitude in degrees" ), &lat );
00050     opts.addOpt< double >( string( "longitude,w" ), string( "east longitude in degrees" ), &lon );
00051 
00052     bool spherical = false;
00053     opts.addOpt< void >( "spherical,s", "use spherical coords input (default false) ", &spherical );
00054 
00055     double distance = 0.01;
00056     opts.addOpt< double >( string( "distance,d" ), string( "distance " ), &distance );
00057 
00058     opts.parseCommandLine( argc, argv );
00059 
00060     int rank         = 0;
00061     int numProcesses = 1;
00062     std::string readopts;
00063 
00064 #ifdef MOAB_HAVE_MPI
00065     if( MPI_Init( &argc, &argv ) ) return 1;
00066     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00067     MPI_Comm_size( MPI_COMM_WORLD, &numProcesses );
00068 
00069     readopts = string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" );  // we do not have to
00070                                                                              // resolve shared ents
00071 #endif
00072     // Instantiate
00073     Core mb;
00074 
00075     // Load the file into a new file set
00076     EntityHandle fileSet;
00077     ErrorCode rval = mb.create_meshset( MESHSET_SET, fileSet );MB_CHK_SET_ERR( rval, "Error creating file set" );
00078     rval = mb.load_file( inputFile.c_str(), &fileSet, readopts.c_str() );MB_CHK_SET_ERR( rval, "Error loading file" );
00079 
00080     if( !rank )
00081     {
00082         cout << " reading file " << inputFile << " on " << numProcesses << " task";
00083         if( numProcesses > 1 ) cout << "s";
00084         cout << "\n";
00085     }
00086     // Get all 2d elements in the file set
00087     Range elems;
00088     rval = mb.get_entities_by_dimension( fileSet, 2, elems );MB_CHK_SET_ERR( rval, "Error getting 2d elements" );
00089 
00090     // create a meshset with close elements
00091     EntityHandle outSet;
00092 
00093     // create meshset
00094     rval = mb.create_meshset( MESHSET_SET, outSet );MB_CHK_ERR( rval );
00095 
00096     // double sphere radius is 1
00097     CartVect point( x, y, z );
00098     if( spherical )
00099     {
00100         std::cout << " use spherical coordinates on the sphere of radius 1.\n";
00101         IntxUtils::SphereCoords pos;
00102         pos.R   = 1;                 // always on the
00103         pos.lat = lat * M_PI / 180;  // defined in math.h
00104         pos.lon = lon * M_PI / 180;
00105         point   = IntxUtils::spherical_to_cart( pos );
00106     }
00107 
00108     Range closeByCells;
00109     for( Range::iterator it = elems.begin(); it != elems.end(); it++ )
00110     {
00111         EntityHandle cell = *it;
00112         CartVect center;
00113         rval = mb.get_coords( &cell, 1, &( center[0] ) );MB_CHK_SET_ERR( rval, "Can't get cell center coords" );
00114         double dist = ( center - point ).length();
00115         if( dist <= distance ) { closeByCells.insert( cell ); }
00116     }
00117 
00118     rval = mb.add_entities( outSet, closeByCells );MB_CHK_SET_ERR( rval, "Can't add to entity set" );
00119 
00120     int numCells = (int)closeByCells.size();
00121 #ifdef MOAB_HAVE_MPI
00122     // Reduce all of the local sums into the global sum
00123     int globalCells;
00124     MPI_Reduce( &numCells, &globalCells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD );
00125     numCells = globalCells;
00126 #endif
00127     if( numCells == 0 )
00128     {
00129         if( !rank )
00130             std::cout << " no close cells to the point " << point << " at distance less than " << distance << "\n";
00131     }
00132     else
00133     {
00134         if( !rank )
00135             std::cout << " write file " << outFile << " with cells closer than " << distance << " from " << point
00136                       << "\n";
00137         string writeOpts;
00138         if( numProcesses > 1 ) writeOpts = string( "PARALLEL=WRITE_PART;" );
00139         rval = mb.write_file( outFile.c_str(), 0, writeOpts.c_str(), &outSet, 1 );MB_CHK_SET_ERR( rval, "Can't write file" );
00140     }
00141 
00142 #ifdef MOAB_HAVE_MPI
00143     if( MPI_Finalize() ) return 1;
00144 #endif
00145     return 0;
00146 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines