 |
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 ContinentsOnGlobe
2 * Description: read a mesh on a sphere and boundaries of continents and major islands,
3 * and write a boundaries mesh file (bound.vtk) and a file with sets for major continents
4 * (map.h5m). Boundaries exist as 2 files, a list of boundary points and a list of loops,
5 * representing each continent and major islands; there are 796 loops (islands). The first one has
6 * 1239 edges/segments (Asia+Europe), while the last one has only 3. It must be a small island :)
7 * ContinentsOnGlobe <input.h5m>
8 * default values: poly2000.h5m : a mesh with 2000 polygons on a sphere of radius 1
9 */
10
11 #include "moab/Core.hpp"
12 #include "moab/ReadUtilIface.hpp"
13 #include <cmath>
14 #include "moab/CartVect.hpp"
15
16 #include <iostream>
17 #include <fstream>
18
19 using namespace moab;
20 using namespace std;
21
22 string bd_name = string( MESH_DIR ) + string( "/../examples/earth/boundary_points.dat" );
23 string loops = string( MESH_DIR ) + string( "/../examples/earth/SaveLoopCounts" );
24 string input_file = string( MESH_DIR ) + string( "/../examples/earth/poly2000.h5m" );
25
26 double getLat( CartVect p )
27 {
28 p.normalize();
29 return asin( p[2] );
30 }
31 double getLon( CartVect p )
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 }
47
48 // we project sphere points on a 2d map. We decide if a point is in interior of a loop
49 // with the span angle test; it will work for any island except Antarctica
50
51 bool interior_point( vector< double >& coords, int& startLoop, int& endLoop, double lat, double lon )
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 }
72
73 int main( int argc, char** argv )
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
148 CartVect center;
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 }
1.9.1.