|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#include "moab/Core.hpp"#include "moab/ReadUtilIface.hpp"#include <cmath>#include "moab/CartVect.hpp"#include <iostream>#include <fstream>
Include dependency graph for ContinentsOnGlobe.cpp:Go to the source code of this file.
Functions | |
| double | getLat (CartVect p) |
| double | getLon (CartVect p) |
| bool | interior_point (vector< double > &coords, int &startLoop, int &endLoop, double lat, double lon) |
| int | main (int argc, char **argv) |
Variables | |
| string | bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" ) |
| string | loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" ) |
| string | input_file = string( MESH_DIR ) + string( "/../examples/earth/poly2000.h5m" ) |
Definition at line 26 of file ContinentsOnGlobe.cpp.
References moab::CartVect::normalize().
{
p.normalize();
return asin( p[2] );
}
Definition at line 31 of file ContinentsOnGlobe.cpp.
References moab::CartVect::normalize().
{
p.normalize();
double lon;
lon = atan2( p[1], p[0] );
if( lon < -2.95 )
{ // separate Asia/Europe from Americas
return 2.0 * M_PI + lon;
}
else
{
return lon;
}
}
| bool interior_point | ( | vector< double > & | coords, |
| int & | startLoop, | ||
| int & | endLoop, | ||
| double | lat, | ||
| double | lon | ||
| ) |
Definition at line 51 of file ContinentsOnGlobe.cpp.
References moab::angle(), and moab::cross().
Referenced by main().
{
// compute oriented angle between point and edges in the loop
// if angle ~ +/-2Pi, it is in the interior
double totalAngle = 0;
for( int i = startLoop; i <= endLoop; i++ )
{
double x1 = coords[2 * i], y1 = coords[2 * i + 1];
int i2 = i + 1;
if( i == endLoop ) i2 = startLoop;
double x2 = coords[2 * i2], y2 = coords[2 * i2 + 1];
CartVect P1( x1 - lat, y1 - lon, 0. );
CartVect P2( x2 - lat, y2 - lon, 0. );
CartVect cross = P1 * P2;
double angle1 = angle( P1, P2 );
if( cross[2] < 0 ) angle1 = -angle1;
totalAngle += angle1;
}
if( fabs( totalAngle + 2 * M_PI ) < 0.1 || fabs( totalAngle - 2 * M_PI ) < 0.1 ) return true;
return false;
}
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
get loops for boundaries
Definition at line 73 of file ContinentsOnGlobe.cpp.
References moab::Interface::add_entities(), bd_name, moab::Range::begin(), center(), moab::Interface::create_meshset(), moab::Range::end(), ErrorCode, moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), getLat(), getLon(), input_file, interior_point(), moab::Interface::load_file(), loops, mb, MB_CHK_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_INTEGER, MESHSET_SET, moab::Range::size(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_data(), and moab::Interface::write_file().
{
// Get MOAB instance
Interface* mb = new( std::nothrow ) Core;
if( NULL == mb ) return 1;
if( argc > 1 )
{
// User has input file for mesh file
input_file = argv[1];
}
cout << "input file : " << input_file << "\n";
std::vector< double > coords;
ifstream bdFile( bd_name.c_str() );
if( !bdFile )
{
cout << endl << "Failed to open file " << bd_name << endl;
;
return 1;
}
coords.resize( 0 );
while( !bdFile.eof() )
{
double co;
bdFile >> co;
coords.push_back( co );
}
cout << " number of boundary points:" << coords.size() / 3 << "\n";
/// get loops for boundaries
vector< int > loopsindx;
ifstream loopFile( loops.c_str() );
while( !loopFile.eof() )
{
int indx;
loopFile >> indx;
loopsindx.push_back( indx );
}
cout << "number of loops: " << loopsindx.size() / 2 << "\n";
// convert loops to lat-lon loops
// it will be easier to determine in-out relations for points in 2d plane
// careful with the international time zone; ( + - Pi longitude )
vector< double > coords2d;
coords2d.resize( coords.size() / 3 * 2 );
for( int i = 0; i < (int)coords.size() / 3; i++ )
{
CartVect p( coords[3 * i], coords[3 * i + 1], coords[3 * i + 2] );
coords2d[2 * i] = getLat( p );
coords2d[2 * i + 1] = getLon( p );
}
ErrorCode rval = mb->load_file( input_file.c_str() );MB_CHK_SET_ERR( rval, "Can't load file" );
// look at the center of element, and see if it is inside the loop
Range cells;
rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
cout << "number of cells: " << cells.size() << "\n";
// tag for continents
Tag tag1;
int defa = -1;
rval = mb->tag_get_handle( "continent", 1, MB_TYPE_INTEGER, tag1, MB_TAG_DENSE | MB_TAG_CREAT, &defa );MB_CHK_SET_ERR( rval, "Trouble creating continent tag" );
EntityHandle islandSets[6];
for( int loop_index = 0; loop_index < 6; loop_index++ )
{
rval = mb->create_meshset( MESHSET_SET, islandSets[loop_index] );MB_CHK_SET_ERR( rval, "Can't create island set" );
int startLoop = loopsindx[2 * loop_index];
int endLoop = loopsindx[2 * loop_index + 1];
vector< EntityHandle > interiorCells;
for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
{
EntityHandle cell = *cit;
// see if it is in the interior of the loop
CartVect center;
rval = mb->get_coords( &cell, 1, &( center[0] ) );MB_CHK_SET_ERR( rval, "Can't get cell center coords" );
double lat = getLat( center ), lon = getLon( center );
// if NA, use some boxes too, for lat/lon
// if (NorthAmerica && (lat < 0.15 || lon < M_PI || lon > 2*M_PI - 0.5) )
// continue;
if( interior_point( coords2d, startLoop, endLoop, lat, lon ) )
{
interiorCells.push_back( cell );
rval = mb->tag_set_data( tag1, &cell, 1, &loop_index );MB_CHK_SET_ERR( rval, "Can't get tag on cell" );
}
}
rval = mb->add_entities( islandSets[loop_index], &interiorCells[0], interiorCells.size() );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
}
std::stringstream islandFile;
islandFile << "map.h5m";
rval = mb->write_file( islandFile.str().c_str() );MB_CHK_SET_ERR( rval, "Can't write island file" );
return 0;
}
Definition at line 22 of file ContinentsOnGlobe.cpp.
| string input_file = string( MESH_DIR ) + string( "/../examples/earth/poly2000.h5m" ) |
Definition at line 24 of file ContinentsOnGlobe.cpp.
Definition at line 23 of file ContinentsOnGlobe.cpp.
1.7.6.1