 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 }
1.9.1.