#include <iostream>
#include <cstdlib>
#include <cstdio>
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CartVect.hpp"
#include "moab/IntxMesh/IntxUtils.hpp"
Go to the source code of this file.
Functions | |
int | main (int argc, char **argv) |
int main | ( | int | argc, |
char ** | argv | ||
) |
Definition at line 30 of file ExtractClose.cpp.
31 {
32 ProgOptions opts;
33
34 // vertex 25 in filt out (-0.784768, 0.594952, -0.173698)
35 double x = -0.784768, y = 0.594952, z = -0.173698;
36 double lat = 28.6551, lon = 61.0204; // from Charlie Zender's page
37 // https://acme-climate.atlassian.net/wiki/spaces/ED/pages/1347092547/Assessing+and+Addressing+Regridding+Weights+in+Map-Files
38
39 string inputFile = test_file_name;
40 opts.addOpt< string >( "inFile,i", "Specify the input file name string ", &inputFile );
41
42 string outFile = "out.h5m";
43 opts.addOpt< string >( "outFile,o", "Specify the output file name string ", &outFile );
44
45 opts.addOpt< double >( string( "xpos,x" ), string( "x position" ), &x );
46 opts.addOpt< double >( string( "ypos,y" ), string( "y position" ), &y );
47 opts.addOpt< double >( string( "zpos,z" ), string( "z position" ), &z );
48
49 opts.addOpt< double >( string( "latitude,t" ), string( "north latitude in degrees" ), &lat );
50 opts.addOpt< double >( string( "longitude,w" ), string( "east longitude in degrees" ), &lon );
51
52 bool spherical = false;
53 opts.addOpt< void >( "spherical,s", "use spherical coords input (default false) ", &spherical );
54
55 double distance = 0.01;
56 opts.addOpt< double >( string( "distance,d" ), string( "distance " ), &distance );
57
58 opts.parseCommandLine( argc, argv );
59
60 int rank = 0;
61 int numProcesses = 1;
62 std::string readopts;
63
64 #ifdef MOAB_HAVE_MPI
65 if( MPI_Init( &argc, &argv ) ) return 1;
66 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
67 MPI_Comm_size( MPI_COMM_WORLD, &numProcesses );
68
69 readopts = string( "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION" ); // we do not have to
70 // resolve shared ents
71 #endif
72 // Instantiate
73 Core mb;
74
75 // Load the file into a new file set
76 EntityHandle fileSet;
77 ErrorCode rval = mb.create_meshset( MESHSET_SET, fileSet );MB_CHK_SET_ERR( rval, "Error creating file set" );
78 rval = mb.load_file( inputFile.c_str(), &fileSet, readopts.c_str() );MB_CHK_SET_ERR( rval, "Error loading file" );
79
80 if( !rank )
81 {
82 cout << " reading file " << inputFile << " on " << numProcesses << " task";
83 if( numProcesses > 1 ) cout << "s";
84 cout << "\n";
85 }
86 // Get all 2d elements in the file set
87 Range elems;
88 rval = mb.get_entities_by_dimension( fileSet, 2, elems );MB_CHK_SET_ERR( rval, "Error getting 2d elements" );
89
90 // create a meshset with close elements
91 EntityHandle outSet;
92
93 // create meshset
94 rval = mb.create_meshset( MESHSET_SET, outSet );MB_CHK_ERR( rval );
95
96 // double sphere radius is 1
97 CartVect point( x, y, z );
98 if( spherical )
99 {
100 std::cout << " use spherical coordinates on the sphere of radius 1.\n";
101 IntxUtils::SphereCoords pos;
102 pos.R = 1; // always on the
103 pos.lat = lat * M_PI / 180; // defined in math.h
104 pos.lon = lon * M_PI / 180;
105 point = IntxUtils::spherical_to_cart( pos );
106 }
107
108 Range closeByCells;
109 for( Range::iterator it = elems.begin(); it != elems.end(); it++ )
110 {
111 EntityHandle cell = *it;
112 CartVect center;
113 rval = mb.get_coords( &cell, 1, &( center[0] ) );MB_CHK_SET_ERR( rval, "Can't get cell center coords" );
114 double dist = ( center - point ).length();
115 if( dist <= distance )
116 {
117 closeByCells.insert( cell );
118 }
119 }
120
121 rval = mb.add_entities( outSet, closeByCells );MB_CHK_SET_ERR( rval, "Can't add to entity set" );
122
123 int numCells = (int)closeByCells.size();
124 #ifdef MOAB_HAVE_MPI
125 // Reduce all of the local sums into the global sum
126 int globalCells;
127 MPI_Reduce( &numCells, &globalCells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD );
128 numCells = globalCells;
129 #endif
130 if( numCells == 0 )
131 {
132 if( !rank )
133 std::cout << " no close cells to the point " << point << " at distance less than " << distance << "\n";
134 }
135 else
136 {
137 if( !rank )
138 std::cout << " write file " << outFile << " with cells closer than " << distance << " from " << point
139 << "\n";
140 string writeOpts;
141 if( numProcesses > 1 ) writeOpts = string( "PARALLEL=WRITE_PART;" );
142 rval = mb.write_file( outFile.c_str(), 0, writeOpts.c_str(), &outSet, 1 );MB_CHK_SET_ERR( rval, "Can't write file" );
143 }
144
145 #ifdef MOAB_HAVE_MPI
146 if( MPI_Finalize() ) return 1;
147 #endif
148 return 0;
149 }
References moab::Core::add_entities(), ProgOptions::addOpt(), moab::Range::begin(), center(), moab::Core::create_meshset(), moab::Range::end(), ErrorCode, moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::Range::insert(), moab::IntxUtils::SphereCoords::lat, length(), moab::Core::load_file(), moab::IntxUtils::SphereCoords::lon, mb, MB_CHK_ERR, MB_CHK_SET_ERR, MESHSET_SET, ProgOptions::parseCommandLine(), moab::IntxUtils::SphereCoords::R, moab::Range::size(), moab::IntxUtils::spherical_to_cart(), test_file_name, and moab::Core::write_file().