Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
BoundaryContinents.cpp
Go to the documentation of this file.
1 /** @example BoundaryContinents
2  * Description: read boundary points and loops that form islands and continents
3  and create an edge mesh file \n
4  *
5  * BoundaryContinents <boundary_points.dat> <SaveLoopCounts>
6  * (default values can run if users don't specify input files)
7  */
8 
9 #include "moab/Core.hpp"
10 #include "moab/ReadUtilIface.hpp"
11 #include "moab/CartVect.hpp"
12 #include <iostream>
13 #include <fstream>
14 
15 using namespace moab;
16 using namespace std;
17 
18 string bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" );
19 string loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" );
20 
21 double getLat( CartVect p )
22 {
23  p.normalize();
24  return asin( p[2] );
25 }
26 double getLon( CartVect p )
27 {
28  p.normalize();
29  double lon;
30 
31  lon = atan2( p[1], p[0] );
32  if( lon < -2.95 )
33  { // this is Behring Strait :)
34  return 2.0 * M_PI + lon;
35  }
36  else
37  {
38  return lon;
39  }
40 }
41 
42 int main( int argc, char** argv )
43 {
44  // Get MOAB instance
45  Interface* mb = new( std::nothrow ) Core;
46  if( NULL == mb ) return 1;
47 
48  if( argc > 1 )
49  {
50  // User has input file for boundary describing continents
51  bd_name = argv[1];
52  }
53  if( argc > 2 )
54  {
55  // User has input file for loops for continents/islands
56  loops = argv[2];
57  }
58 
59  std::vector< double > coords;
60  ifstream bdFile( bd_name.c_str() );
61  if( !bdFile )
62  {
63  cout << endl << "Failed to open file " << bd_name << endl;
64  ;
65  return 1;
66  }
67  coords.resize( 0 );
68  while( !bdFile.eof() )
69  {
70  double co;
71  bdFile >> co;
72  coords.push_back( co );
73  }
74  cout << " number of boundary points:" << coords.size() / 3 << "\n";
75  /// get loops for boundaries
76  vector< int > loopsindx;
77  ifstream loopFile( loops.c_str() );
78  while( !loopFile.eof() )
79  {
80  int indx;
81  loopFile >> indx;
82  loopsindx.push_back( indx );
83  }
84  cout << "number of loops: " << loopsindx.size() / 2 << "\n";
85  ReadUtilIface* readMeshIface;
86  mb->query_interface( readMeshIface );
87  Range verts;
88 
89  ErrorCode rval = mb->create_vertices( &coords[0], coords.size() / 3, verts );MB_CHK_SET_ERR( rval, "do not create boundary vertices" );
90 
91  Tag gid = mb->globalId_tag();
92 
93  int num_edges = 0;
94  for( size_t i = 0; i < loopsindx.size() / 2; i++ )
95  {
96  num_edges += ( loopsindx[2 * i + 1] - loopsindx[2 * i] + 1 );
97  }
98  EntityHandle handle;
99  EntityHandle* conn_array;
100  rval = readMeshIface->get_element_connect( num_edges, 2, MBEDGE, 1, handle, conn_array );MB_CHK_SET_ERR( rval, "do not elems" );
101  int i1 = 0; // index in edge connectivity
102  for( size_t i = 0; i < loopsindx.size() / 2; i++ )
103  {
104  for( int j = loopsindx[2 * i]; j <= loopsindx[2 * i + 1]; j++ )
105  {
106  int j2;
107  j2 = j;
108  if( j == loopsindx[2 * i + 1] ) j2 = loopsindx[2 * i] - 1; // close the loop
109  conn_array[2 * i1] = verts[j - 1]; // vertices for sure start at 1
110  conn_array[2 * i1 + 1] = verts[j2]; // vertices for sure start at 1
111  i1++;
112  }
113  // the last vertex is first one in the loop
114  }
115 
116  Range bedges( handle, handle + num_edges - 1 );
117  EntityHandle bound_set;
118  rval = mb->create_meshset( MESHSET_SET, bound_set );MB_CHK_SET_ERR( rval, "Can't create boundary edges set" );
119 
120  rval = mb->add_entities( bound_set, bedges );MB_CHK_SET_ERR( rval, "Can't add edges to boundary set" );
121 
122  // set global ids for vertices and edges
123  vector< int > gids;
124  gids.resize( verts.size() );
125  for( int j = 0; j < (int)verts.size(); j++ )
126  {
127  gids[j] = j + 1;
128  }
129  rval = mb->tag_set_data( gid, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids on verts" );
130  // we have less edges than vertices, we can reuse the array for edges too
131  rval = mb->tag_set_data( gid, bedges, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global ids on edges" );
132 
133  mb->write_file( "bound.vtk", 0, 0, &bound_set, 1 );
134 
135  // now project the mesh in 2d plane
136  // get all the vertices coordinates
137  std::vector< CartVect > co3;
138  co3.resize( verts.size() );
139  rval = mb->get_coords( verts, &( co3[0][0] ) );MB_CHK_SET_ERR( rval, "Can't get vertex coords" );
140  for( size_t i = 0; i < verts.size(); i++ )
141  {
142  CartVect p = co3[i];
143  double lat1 = getLat( p );
144  double lon1 = getLon( p );
145  co3[i] = CartVect( lon1, lat1, 0. );
146  }
147  rval = mb->set_coords( verts, &( co3[0][0] ) );MB_CHK_SET_ERR( rval, "Can't set new vertex coords" );
148  // remove edges in 2d that are too long (longer than 6; they are on the other side..., periodic)
149 
150  Range longEdges;
151  for( Range::iterator eit = bedges.begin(); eit != bedges.end(); ++eit )
152  {
153  EntityHandle eh = *eit;
154  const EntityHandle* conn;
155  int nconn;
156  rval = mb->get_connectivity( eh, conn, nconn );MB_CHK_ERR( rval );
157  // get coordinates of nodes
158  CartVect c2[2];
159  rval = mb->get_coords( conn, 2, &( c2[0][0] ) );MB_CHK_ERR( rval );
160  double dist = ( c2[1] - c2[0] ).length_squared();
161  if( dist > 36. ) // too long edge in 2d, remove it
162  longEdges.insert( eh );
163  }
164  rval = mb->delete_entities( longEdges );MB_CHK_ERR( rval );
165 
166  mb->write_file( "bound2d.vtk" );
167  delete mb;
168 
169  return 0;
170 }