MOAB: Mesh Oriented datABase  (version 5.2.1)
BoundaryContinents.cpp
Go to the documentation of this file.
00001 /** @example BoundaryContinents
00002  * Description: read boundary points and loops that form islands and continents
00003       and create an edge mesh file \n
00004  *
00005  *    BoundaryContinents  <boundary_points.dat> <SaveLoopCounts>
00006  * (default values can run if users don't specify input files)
00007  */
00008 
00009 #include "moab/Core.hpp"
00010 #include "moab/ReadUtilIface.hpp"
00011 #include "moab/CartVect.hpp"
00012 #include <iostream>
00013 #include <fstream>
00014 
00015 using namespace moab;
00016 using namespace std;
00017 
00018 string bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" );
00019 string loops   = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" );
00020 
00021 double getLat( CartVect p )
00022 {
00023     p.normalize();
00024     return asin( p[2] );
00025 }
00026 double getLon( CartVect p )
00027 {
00028     p.normalize();
00029     double lon;
00030 
00031     lon = atan2( p[1], p[0] );
00032     if( lon < -2.95 )
00033     {  // this is Behring Strait :)
00034         return 2.0 * M_PI + lon;
00035     }
00036     else
00037     {
00038         return lon;
00039     }
00040 }
00041 
00042 int main( int argc, char** argv )
00043 {
00044     // Get MOAB instance
00045     Interface* mb = new( std::nothrow ) Core;
00046     if( NULL == mb ) return 1;
00047 
00048     if( argc > 1 )
00049     {
00050         // User has input file for boundary describing continents
00051         bd_name = argv[1];
00052     }
00053     if( argc > 2 )
00054     {
00055         // User has input file for loops for continents/islands
00056         loops = argv[2];
00057     }
00058 
00059     std::vector< double > coords;
00060     ifstream bdFile( bd_name.c_str() );
00061     if( !bdFile )
00062     {
00063         cout << endl << "Failed to open file " << bd_name << endl;
00064         ;
00065         return 1;
00066     }
00067     coords.resize( 0 );
00068     while( !bdFile.eof() )
00069     {
00070         double co;
00071         bdFile >> co;
00072         coords.push_back( co );
00073     }
00074     cout << " number of boundary points:" << coords.size() / 3 << "\n";
00075     /// get loops for boundaries
00076     vector< int > loopsindx;
00077     ifstream loopFile( loops.c_str() );
00078     while( !loopFile.eof() )
00079     {
00080         int indx;
00081         loopFile >> indx;
00082         loopsindx.push_back( indx );
00083     }
00084     cout << "number of loops: " << loopsindx.size() / 2 << "\n";
00085     ReadUtilIface* readMeshIface;
00086     mb->query_interface( readMeshIface );
00087     Range verts;
00088 
00089     ErrorCode rval = mb->create_vertices( &coords[0], coords.size() / 3, verts );MB_CHK_SET_ERR( rval, "do not create boundary vertices" );
00090 
00091     Tag gid = mb->globalId_tag();
00092 
00093     int num_edges = 0;
00094     for( size_t i = 0; i < loopsindx.size() / 2; i++ )
00095     {
00096         num_edges += ( loopsindx[2 * i + 1] - loopsindx[2 * i] + 1 );
00097     }
00098     EntityHandle handle;
00099     EntityHandle* conn_array;
00100     rval = readMeshIface->get_element_connect( num_edges, 2, MBEDGE, 1, handle, conn_array );MB_CHK_SET_ERR( rval, "do not elems" );
00101     int i1 = 0;  // index in edge connectivity
00102     for( size_t i = 0; i < loopsindx.size() / 2; i++ )
00103     {
00104         for( int j = loopsindx[2 * i]; j <= loopsindx[2 * i + 1]; j++ )
00105         {
00106             int j2;
00107             j2 = j;
00108             if( j == loopsindx[2 * i + 1] ) j2 = loopsindx[2 * i] - 1;  // close the loop
00109             conn_array[2 * i1]     = verts[j - 1];                      // vertices for sure start at 1
00110             conn_array[2 * i1 + 1] = verts[j2];                         // vertices for sure start at 1
00111             i1++;
00112         }
00113         // the last vertex is first one in the loop
00114     }
00115 
00116     Range bedges( handle, handle + num_edges - 1 );
00117     EntityHandle bound_set;
00118     rval = mb->create_meshset( MESHSET_SET, bound_set );MB_CHK_SET_ERR( rval, "Can't create boundary edges set" );
00119 
00120     rval = mb->add_entities( bound_set, bedges );MB_CHK_SET_ERR( rval, "Can't add edges to boundary set" );
00121 
00122     // set global ids for vertices and edges
00123     vector< int > gids;
00124     gids.resize( verts.size() );
00125     for( int j = 0; j < (int)verts.size(); j++ )
00126     {
00127         gids[j] = j + 1;
00128     }
00129     rval = mb->tag_set_data( gid, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids on verts" );
00130     // we have less edges than vertices, we can reuse the array for edges too
00131     rval = mb->tag_set_data( gid, bedges, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids on edges" );
00132 
00133     mb->write_file( "bound.vtk", 0, 0, &bound_set, 1 );
00134 
00135     // now project the mesh in 2d plane
00136     // get all the vertices coordinates
00137     std::vector< CartVect > co3;
00138     co3.resize( verts.size() );
00139     rval = mb->get_coords( verts, &( co3[0][0] ) );MB_CHK_SET_ERR( rval, "Can't get vertex coords" );
00140     for( size_t i = 0; i < verts.size(); i++ )
00141     {
00142         CartVect p  = co3[i];
00143         double lat1 = getLat( p );
00144         double lon1 = getLon( p );
00145         co3[i]      = CartVect( lon1, lat1, 0. );
00146     }
00147     rval = mb->set_coords( verts, &( co3[0][0] ) );MB_CHK_SET_ERR( rval, "Can't set new vertex coords" );
00148     // remove edges in 2d that are too long (longer than 6; they are on the other side..., periodic)
00149 
00150     Range longEdges;
00151     for( Range::iterator eit = bedges.begin(); eit != bedges.end(); ++eit )
00152     {
00153         EntityHandle eh = *eit;
00154         const EntityHandle* conn;
00155         int nconn;
00156         rval = mb->get_connectivity( eh, conn, nconn );MB_CHK_ERR( rval );
00157         // get coordinates of nodes
00158         CartVect c2[2];
00159         rval = mb->get_coords( conn, 2, &( c2[0][0] ) );MB_CHK_ERR( rval );
00160         double dist = ( c2[1] - c2[0] ).length_squared();
00161         if( dist > 36. )  // too long edge in 2d, remove it
00162             longEdges.insert( eh );
00163     }
00164     rval = mb->delete_entities( longEdges );MB_CHK_ERR( rval );
00165 
00166     mb->write_file( "bound2d.vtk" );
00167     delete mb;
00168 
00169     return 0;
00170 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines