MOAB: Mesh Oriented datABase  (version 5.2.1)
ContinentsOnGlobe.cpp
Go to the documentation of this file.
00001 /** @example ContinentsOnGlobe
00002  * Description: read a mesh on a sphere and boundaries of continents and major islands,
00003  *   and write a boundaries mesh file (bound.vtk) and a file with sets for major continents
00004  * (map.h5m). Boundaries exist as 2 files, a list of boundary points and a list of loops,
00005  * representing each continent and major islands; there are 796 loops (islands). The first one has
00006  * 1239 edges/segments (Asia+Europe), while the last one has only 3. It must be a small island :)
00007  *    ContinentsOnGlobe  <input.h5m>
00008  *  default values: poly2000.h5m : a mesh with 2000 polygons on a sphere of radius 1
00009  */
00010 
00011 #include "moab/Core.hpp"
00012 #include "moab/ReadUtilIface.hpp"
00013 #include "math.h"
00014 #include "moab/CartVect.hpp"
00015 
00016 #include <iostream>
00017 #include <fstream>
00018 
00019 using namespace moab;
00020 using namespace std;
00021 
00022 string bd_name    = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" );
00023 string loops      = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" );
00024 string input_file = string( MESH_DIR ) + string( "/../examples/earth/poly2000.h5m" );
00025 
00026 double getLat( CartVect p )
00027 {
00028     p.normalize();
00029     return asin( p[2] );
00030 }
00031 double getLon( CartVect p )
00032 {
00033     p.normalize();
00034     double lon;
00035 
00036     lon = atan2( p[1], p[0] );
00037 
00038     if( lon < -2.95 )
00039     {  // separate Asia/Europe from Americas
00040         return 2.0 * M_PI + lon;
00041     }
00042     else
00043     {
00044         return lon;
00045     }
00046 }
00047 
00048 // we project sphere points on a 2d map. We decide if a point is in interior of a loop
00049 // with the span angle test; it will work for any island except Antarctica
00050 
00051 bool interior_point( vector< double >& coords, int& startLoop, int& endLoop, double lat, double lon )
00052 {
00053     // compute oriented angle between point and edges in the loop
00054     // if angle ~ +/-2Pi, it is in the interior
00055     double totalAngle = 0;
00056     for( int i = startLoop; i <= endLoop; i++ )
00057     {
00058         double x1 = coords[2 * i], y1 = coords[2 * i + 1];
00059         int i2 = i + 1;
00060         if( i == endLoop ) i2 = startLoop;
00061         double x2 = coords[2 * i2], y2 = coords[2 * i2 + 1];
00062         CartVect P1( x1 - lat, y1 - lon, 0. );
00063         CartVect P2( x2 - lat, y2 - lon, 0. );
00064         CartVect cross = P1 * P2;
00065         double angle1  = angle( P1, P2 );
00066         if( cross[2] < 0 ) angle1 = -angle1;
00067         totalAngle += angle1;
00068     }
00069     if( fabs( totalAngle + 2 * M_PI ) < 0.1 || fabs( totalAngle - 2 * M_PI ) < 0.1 ) return true;
00070     return false;
00071 }
00072 
00073 int main( int argc, char** argv )
00074 {
00075     // Get MOAB instance
00076     Interface* mb = new( std::nothrow ) Core;
00077     if( NULL == mb ) return 1;
00078 
00079     if( argc > 1 )
00080     {
00081         // User has input file for mesh file
00082         input_file = argv[1];
00083     }
00084 
00085     cout << "input file : " << input_file << "\n";
00086     std::vector< double > coords;
00087     ifstream bdFile( bd_name.c_str() );
00088     if( !bdFile )
00089     {
00090         cout << endl << "Failed to open file " << bd_name << endl;
00091         ;
00092         return 1;
00093     }
00094     coords.resize( 0 );
00095     while( !bdFile.eof() )
00096     {
00097         double co;
00098         bdFile >> co;
00099         coords.push_back( co );
00100     }
00101     cout << " number of boundary points:" << coords.size() / 3 << "\n";
00102     /// get loops for boundaries
00103     vector< int > loopsindx;
00104     ifstream loopFile( loops.c_str() );
00105     while( !loopFile.eof() )
00106     {
00107         int indx;
00108         loopFile >> indx;
00109         loopsindx.push_back( indx );
00110     }
00111     cout << "number of loops: " << loopsindx.size() / 2 << "\n";
00112     // convert loops to lat-lon loops
00113     // it will be easier to determine in-out relations for points in 2d plane
00114     // careful with the international time zone; ( + - Pi longitude )
00115     vector< double > coords2d;
00116     coords2d.resize( coords.size() / 3 * 2 );
00117     for( int i = 0; i < (int)coords.size() / 3; i++ )
00118     {
00119         CartVect p( coords[3 * i], coords[3 * i + 1], coords[3 * i + 2] );
00120         coords2d[2 * i]     = getLat( p );
00121         coords2d[2 * i + 1] = getLon( p );
00122     }
00123 
00124     ErrorCode rval = mb->load_file( input_file.c_str() );MB_CHK_SET_ERR( rval, "Can't load file" );
00125     // look at the center of element, and see if it is inside the loop
00126 
00127     Range cells;
00128     rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
00129 
00130     cout << "number of cells: " << cells.size() << "\n";
00131 
00132     // tag for continents
00133     Tag tag1;
00134     int defa = -1;
00135     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" );
00136     EntityHandle islandSets[6];
00137     for( int loop_index = 0; loop_index < 6; loop_index++ )
00138     {
00139         rval = mb->create_meshset( MESHSET_SET, islandSets[loop_index] );MB_CHK_SET_ERR( rval, "Can't create island set" );
00140         int startLoop = loopsindx[2 * loop_index];
00141         int endLoop   = loopsindx[2 * loop_index + 1];
00142 
00143         vector< EntityHandle > interiorCells;
00144         for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
00145         {
00146             EntityHandle cell = *cit;
00147             // see if it is in the interior of the loop
00148             CartVect center;
00149             rval = mb->get_coords( &cell, 1, &( center[0] ) );MB_CHK_SET_ERR( rval, "Can't get cell center coords" );
00150             double lat = getLat( center ), lon = getLon( center );
00151             // if NA, use some boxes too, for lat/lon
00152             // if (NorthAmerica && (lat < 0.15 || lon < M_PI || lon > 2*M_PI - 0.5) )
00153             //  continue;
00154 
00155             if( interior_point( coords2d, startLoop, endLoop, lat, lon ) )
00156             {
00157                 interiorCells.push_back( cell );
00158                 rval = mb->tag_set_data( tag1, &cell, 1, &loop_index );MB_CHK_SET_ERR( rval, "Can't get tag on cell" );
00159             }
00160         }
00161 
00162         rval = mb->add_entities( islandSets[loop_index], &interiorCells[0], interiorCells.size() );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
00163     }
00164 
00165     std::stringstream islandFile;
00166 
00167     islandFile << "map.h5m";
00168 
00169     rval = mb->write_file( islandFile.str().c_str() );MB_CHK_SET_ERR( rval, "Can't write island file" );
00170     return 0;
00171 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines