MOAB: Mesh Oriented datABase  (version 5.2.1)
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     int fail = MPI_Init( &argc, &argv );
00066     if( fail ) return 1;
00067     MPI_Comm_rank( MPI_COMM_WORLD, &rank );
00068     MPI_Comm_size( MPI_COMM_WORLD, &numProcesses );
00069 
00070     readopts = string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" );  // we do not have to
00071                                                                              // resolve shared ents
00072 #endif
00073     // Instantiate
00074     Core mb;
00075 
00076     // Load the file into a new file set
00077     EntityHandle fileSet;
00078     ErrorCode rval = mb.create_meshset( MESHSET_SET, fileSet );MB_CHK_SET_ERR( rval, "Error creating file set" );
00079     rval = mb.load_file( inputFile.c_str(), &fileSet, readopts.c_str() );MB_CHK_SET_ERR( rval, "Error loading file" );
00080 
00081     if( !rank )
00082     {
00083         cout << " reading file " << inputFile << " on " << numProcesses << " task";
00084         if( numProcesses > 1 ) cout << "s";
00085         cout << "\n";
00086     }
00087     // Get all 2d elements in the file set
00088     Range elems;
00089     rval = mb.get_entities_by_dimension( fileSet, 2, elems );MB_CHK_SET_ERR( rval, "Error getting 2d elements" );
00090 
00091     // create a meshset with close elements
00092     EntityHandle outSet;
00093 
00094     // create meshset
00095     rval = mb.create_meshset( MESHSET_SET, outSet );MB_CHK_ERR( rval );
00096 
00097     // double sphere radius is 1
00098     CartVect point( x, y, z );
00099     if( spherical )
00100     {
00101         std::cout << " use spherical coordinates on the sphere of radius 1.\n";
00102         IntxUtils::SphereCoords pos;
00103         pos.R   = 1;                 // always on the
00104         pos.lat = lat * M_PI / 180;  // defined in math.h
00105         pos.lon = lon * M_PI / 180;
00106         point   = IntxUtils::spherical_to_cart( pos );
00107     }
00108 
00109     Range closeByCells;
00110     for( Range::iterator it = elems.begin(); it != elems.end(); it++ )
00111     {
00112         EntityHandle cell = *it;
00113         CartVect center;
00114         rval = mb.get_coords( &cell, 1, &( center[0] ) );MB_CHK_SET_ERR( rval, "Can't get cell center coords" );
00115         double dist = ( center - point ).length();
00116         if( dist <= distance ) { closeByCells.insert( cell ); }
00117     }
00118 
00119     rval = mb.add_entities( outSet, closeByCells );MB_CHK_SET_ERR( rval, "Can't add to entity set" );
00120 
00121     int numCells = (int)closeByCells.size();
00122 #ifdef MOAB_HAVE_MPI
00123     // Reduce all of the local sums into the global sum
00124     int globalCells;
00125     MPI_Reduce( &numCells, &globalCells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD );
00126     numCells = globalCells;
00127 #endif
00128     if( numCells == 0 )
00129     {
00130         if( !rank )
00131             std::cout << " no close cells to the point " << point << " at distance less than " << distance << "\n";
00132     }
00133     else
00134     {
00135         if( !rank )
00136             std::cout << " write file " << outFile << " with cells closer than " << distance << " from " << point
00137                       << "\n";
00138         string writeOpts;
00139         if( numProcesses > 1 ) writeOpts = string( "PARALLEL=WRITE_PART;" );
00140         rval = mb.write_file( outFile.c_str(), 0, writeOpts.c_str(), &outSet, 1 );MB_CHK_SET_ERR( rval, "Can't write file" );
00141     }
00142     return 0;
00143 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines