BoundaryContinents.cpp

Go to the documentation of this file.

/** @example BoundaryContinents
 * Description: read boundary points and loops that form islands and continents
 and create an edge mesh file \n
 *
 * BoundaryContinents <boundary_points.dat> <SaveLoopCounts>
 * (default values can run if users don't specify input files)
 */
 
#include "moab/Core.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/CartVect.hpp"
#include <iostream>
#include <fstream>
 
using namespace moab;
using namespace std;
 
string bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" );
string loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" );
 
double getLat( CartVect p )
{
 p.normalize();
 return asin( p[2] );
}
double getLon( CartVect p )
{
 p.normalize();
 double lon;
 
 lon = atan2( p[1], p[0] );
 if( lon < -2.95 )
 { // this is Behring Strait :)
 return 2.0 * M_PI + lon;
 }
 else
 {
 return lon;
 }
}
 
int main( int argc, char** argv )
{
 // Get MOAB instance
 Interface* mb = new( std::nothrow ) Core;
 if( NULL == mb ) return 1;
 
 if( argc > 1 )
 {
 // User has input file for boundary describing continents
 bd_name = argv[1];
 }
 if( argc > 2 )
 {
 // User has input file for loops for continents/islands
 loops = argv[2];
 }
 
 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";
 ReadUtilIface* readMeshIface;
 mb->query_interface( readMeshIface );
 Range verts;
 
 ErrorCode rval = mb->create_vertices( &coords[0], coords.size() / 3, verts );MB_CHK_SET_ERR( rval, "do not create boundary vertices" );
 
 Tag gid = mb->globalId_tag();
 
 int num_edges = 0;
 for( size_t i = 0; i < loopsindx.size() / 2; i++ )
 {
 num_edges += ( loopsindx[2 * i + 1] - loopsindx[2 * i] + 1 );
 }
 EntityHandle handle;
 EntityHandle* conn_array;
 rval = readMeshIface->get_element_connect( num_edges, 2, MBEDGE, 1, handle, conn_array );MB_CHK_SET_ERR( rval, "do not elems" );
 int i1 = 0; // index in edge connectivity
 for( size_t i = 0; i < loopsindx.size() / 2; i++ )
 {
 for( int j = loopsindx[2 * i]; j <= loopsindx[2 * i + 1]; j++ )
 {
 int j2;
 j2 = j;
 if( j == loopsindx[2 * i + 1] ) j2 = loopsindx[2 * i] - 1; // close the loop
 conn_array[2 * i1] = verts[j - 1]; // vertices for sure start at 1
 conn_array[2 * i1 + 1] = verts[j2]; // vertices for sure start at 1
 i1++;
 }
 // the last vertex is first one in the loop
 }
 
 Range bedges( handle, handle + num_edges - 1 );
 EntityHandle bound_set;
 rval = mb->create_meshset( MESHSET_SET, bound_set );MB_CHK_SET_ERR( rval, "Can't create boundary edges set" );
 
 rval = mb->add_entities( bound_set, bedges );MB_CHK_SET_ERR( rval, "Can't add edges to boundary set" );
 
 // set global ids for vertices and edges
 vector< int > gids;
 gids.resize( verts.size() );
 for( int j = 0; j < (int)verts.size(); j++ )
 {
 gids[j] = j + 1;
 }
 rval = mb->tag_set_data( gid, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids on verts" );
 // we have less edges than vertices, we can reuse the array for edges too
 rval = mb->tag_set_data( gid, bedges, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids on edges" );
 
 mb->write_file( "bound.vtk", 0, 0, &bound_set, 1 );
 
 // now project the mesh in 2d plane
 // get all the vertices coordinates
 std::vector< CartVect > co3;
 co3.resize( verts.size() );
 rval = mb->get_coords( verts, &( co3[0][0] ) );MB_CHK_SET_ERR( rval, "Can't get vertex coords" );
 for( size_t i = 0; i < verts.size(); i++ )
 {
 CartVect p = co3[i];
 double lat1 = getLat( p );
 double lon1 = getLon( p );
 co3[i] = CartVect( lon1, lat1, 0. );
 }
 rval = mb->set_coords( verts, &( co3[0][0] ) );MB_CHK_SET_ERR( rval, "Can't set new vertex coords" );
 // remove edges in 2d that are too long (longer than 6; they are on the other side..., periodic)
 
 Range longEdges;
 for( Range::iterator eit = bedges.begin(); eit != bedges.end(); ++eit )
 {
 EntityHandle eh = *eit;
 const EntityHandle* conn;
 int nconn;
 rval = mb->get_connectivity( eh, conn, nconn );MB_CHK_ERR( rval );
 // get coordinates of nodes
 CartVect c2[2];
 rval = mb->get_coords( conn, 2, &( c2[0][0] ) );MB_CHK_ERR( rval );
 double dist = ( c2[1] - c2[0] ).length_squared();
 if( dist > 36. ) // too long edge in 2d, remove it
 longEdges.insert( eh );
 }
 rval = mb->delete_entities( longEdges );MB_CHK_ERR( rval );
 
 mb->write_file( "bound2d.vtk" );
 delete mb;
 
 return 0;
}