Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ContinentsOnGlobe.cpp File Reference
#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" )
 

Function Documentation

◆ getLat()

double getLat ( CartVect  p)

Definition at line 26 of file ContinentsOnGlobe.cpp.

27 {
28  p.normalize();
29  return asin( p[2] );
30 }

References moab::CartVect::normalize().

Referenced by main().

◆ getLon()

double getLon ( CartVect  p)

Definition at line 31 of file ContinentsOnGlobe.cpp.

32 {
33  p.normalize();
34  double lon;
35 
36  lon = atan2( p[1], p[0] );
37 
38  if( lon < -2.95 )
39  { // separate Asia/Europe from Americas
40  return 2.0 * M_PI + lon;
41  }
42  else
43  {
44  return lon;
45  }
46 }

References moab::CartVect::normalize().

Referenced by main().

◆ interior_point()

bool interior_point ( vector< double > &  coords,
int &  startLoop,
int &  endLoop,
double  lat,
double  lon 
)

Definition at line 51 of file ContinentsOnGlobe.cpp.

52 {
53  // compute oriented angle between point and edges in the loop
54  // if angle ~ +/-2Pi, it is in the interior
55  double totalAngle = 0;
56  for( int i = startLoop; i <= endLoop; i++ )
57  {
58  double x1 = coords[2 * i], y1 = coords[2 * i + 1];
59  int i2 = i + 1;
60  if( i == endLoop ) i2 = startLoop;
61  double x2 = coords[2 * i2], y2 = coords[2 * i2 + 1];
62  CartVect P1( x1 - lat, y1 - lon, 0. );
63  CartVect P2( x2 - lat, y2 - lon, 0. );
64  CartVect cross = P1 * P2;
65  double angle1 = angle( P1, P2 );
66  if( cross[2] < 0 ) angle1 = -angle1;
67  totalAngle += angle1;
68  }
69  if( fabs( totalAngle + 2 * M_PI ) < 0.1 || fabs( totalAngle - 2 * M_PI ) < 0.1 ) return true;
70  return false;
71 }

References moab::angle(), and moab::cross().

Referenced by main().

◆ main()

int main ( int  argc,
char **  argv 
)

get loops for boundaries

Definition at line 73 of file ContinentsOnGlobe.cpp.

74 {
75  // Get MOAB instance
76  Interface* mb = new( std::nothrow ) Core;
77  if( NULL == mb ) return 1;
78 
79  if( argc > 1 )
80  {
81  // User has input file for mesh file
82  input_file = argv[1];
83  }
84 
85  cout << "input file : " << input_file << "\n";
86  std::vector< double > coords;
87  ifstream bdFile( bd_name.c_str() );
88  if( !bdFile )
89  {
90  cout << endl << "Failed to open file " << bd_name << endl;
91  ;
92  return 1;
93  }
94  coords.resize( 0 );
95  while( !bdFile.eof() )
96  {
97  double co;
98  bdFile >> co;
99  coords.push_back( co );
100  }
101  cout << " number of boundary points:" << coords.size() / 3 << "\n";
102  /// get loops for boundaries
103  vector< int > loopsindx;
104  ifstream loopFile( loops.c_str() );
105  while( !loopFile.eof() )
106  {
107  int indx;
108  loopFile >> indx;
109  loopsindx.push_back( indx );
110  }
111  cout << "number of loops: " << loopsindx.size() / 2 << "\n";
112  // convert loops to lat-lon loops
113  // it will be easier to determine in-out relations for points in 2d plane
114  // careful with the international time zone; ( + - Pi longitude )
115  vector< double > coords2d;
116  coords2d.resize( coords.size() / 3 * 2 );
117  for( int i = 0; i < (int)coords.size() / 3; i++ )
118  {
119  CartVect p( coords[3 * i], coords[3 * i + 1], coords[3 * i + 2] );
120  coords2d[2 * i] = getLat( p );
121  coords2d[2 * i + 1] = getLon( p );
122  }
123 
124  ErrorCode rval = mb->load_file( input_file.c_str() );MB_CHK_SET_ERR( rval, "Can't load file" );
125  // look at the center of element, and see if it is inside the loop
126 
127  Range cells;
128  rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "Can't get cells" );
129 
130  cout << "number of cells: " << cells.size() << "\n";
131 
132  // tag for continents
133  Tag tag1;
134  int defa = -1;
135  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" );
136  EntityHandle islandSets[6];
137  for( int loop_index = 0; loop_index < 6; loop_index++ )
138  {
139  rval = mb->create_meshset( MESHSET_SET, islandSets[loop_index] );MB_CHK_SET_ERR( rval, "Can't create island set" );
140  int startLoop = loopsindx[2 * loop_index];
141  int endLoop = loopsindx[2 * loop_index + 1];
142 
143  vector< EntityHandle > interiorCells;
144  for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
145  {
146  EntityHandle cell = *cit;
147  // see if it is in the interior of the loop
149  rval = mb->get_coords( &cell, 1, &( center[0] ) );MB_CHK_SET_ERR( rval, "Can't get cell center coords" );
150  double lat = getLat( center ), lon = getLon( center );
151  // if NA, use some boxes too, for lat/lon
152  // if (NorthAmerica && (lat < 0.15 || lon < M_PI || lon > 2*M_PI - 0.5) )
153  // continue;
154 
155  if( interior_point( coords2d, startLoop, endLoop, lat, lon ) )
156  {
157  interiorCells.push_back( cell );
158  rval = mb->tag_set_data( tag1, &cell, 1, &loop_index );MB_CHK_SET_ERR( rval, "Can't get tag on cell" );
159  }
160  }
161 
162  rval = mb->add_entities( islandSets[loop_index], &interiorCells[0], interiorCells.size() );MB_CHK_SET_ERR( rval, "Can't add entities to set" );
163  }
164 
165  std::stringstream islandFile;
166 
167  islandFile << "map.h5m";
168 
169  rval = mb->write_file( islandFile.str().c_str() );MB_CHK_SET_ERR( rval, "Can't write island file" );
170  return 0;
171 }

References moab::Core::add_entities(), bd_name, moab::Range::begin(), center(), moab::Core::create_meshset(), moab::Range::end(), ErrorCode, moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), getLat(), getLon(), input_file, interior_point(), moab::Core::load_file(), loops, mb, MB_CHK_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_INTEGER, MESHSET_SET, moab::Range::size(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), and moab::Core::write_file().

Variable Documentation

◆ bd_name

string bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" )

Definition at line 22 of file ContinentsOnGlobe.cpp.

Referenced by main().

◆ input_file

◆ loops

string loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" )

Definition at line 23 of file ContinentsOnGlobe.cpp.

Referenced by main().