1 #include "moab/Core.hpp"
2 #include "moab/Range.hpp"
3 #include "MBTagConventions.hpp"
4 #include <iostream>
5 #include <fstream>
6 #include <limits>
7 #include <cstdlib>
8 #include <math.h>
9
10 /* Exit values */
11 #define USAGE_ERROR 1
12 #define READ_ERROR 2
13 #define WRITE_ERROR 3
14 #define SURFACE_NOT_FOUND 4
15 #define OTHER_ERROR 5
16
usage_error(const char * name)17 static void usage_error( const char* name )
18 {
19 std::cerr << "Usage: " << name <<
20 " [-g|-p] <Surface_ID> <input_file>" << std::endl
21 << "\t-g - Write GNU Plot data file (default)." << std::endl
22 << "\t-p - Write encapsulated postscript file" << std::endl
23 << "\t-s - Write an SVG file" << std::endl
24 << "\t<Surface_ID> - ID of surface containing mesh to export (0 for entire file)." << std::endl
25 << "\t<input_file> - Mesh file to read." << std::endl
26 << std::endl
27 << " This utility plots the mesh of a single geometric surface "
28 << "projected to a plane. The output file is written to stdout."
29 << std::endl;
30
31 exit(USAGE_ERROR);
32 }
33
34 struct CartVect3D {
35
36 double x, y, z;
37
CartVect3DCartVect3D38 CartVect3D() : x(0.0), y(0.0), z(0.0) {}
39
CartVect3DCartVect3D40 CartVect3D( double x_, double y_, double z_ )
41 : x(x_), y(y_), z(z_) {}
42
operator +=CartVect3D43 CartVect3D& operator+=( const CartVect3D& o )
44 { x += o.x; y += o.y; z += o.z; return *this; }
45
operator -=CartVect3D46 CartVect3D& operator-=( const CartVect3D& o )
47 { x -= o.x; y -= o.y; z -= o.z; return *this; }
48
49 CartVect3D& operator*=( const CartVect3D& );
50
operator +=CartVect3D51 CartVect3D& operator+=( double v )
52 { x += v; y += v; z += v; return *this; }
53
operator -=CartVect3D54 CartVect3D& operator-=( double v )
55 { x -= v; y -= v; z -= v; return *this; }
56
operator *=CartVect3D57 CartVect3D& operator*=( double v )
58 { x *= v; y *= v; z *= v; return *this; }
59
operator /=CartVect3D60 CartVect3D& operator/=( double v )
61 { x /= v; y /= v; z /= v; return *this; }
62
lenCartVect3D63 double len() const
64 { return sqrt( x*x + y*y + z*z ); }
65 };
66
67 //static CartVect3D operator-( const CartVect3D &a )
68 // { return CartVect3D(-a.z, -a.y, -a.z); }
69
70 //static CartVect3D operator+( const CartVect3D &a, const CartVect3D &b )
71 // { return CartVect3D(a.x+b.x, a.y+b.y, a.z+b.z); }
72
operator -(const CartVect3D & a,const CartVect3D & b)73 static CartVect3D operator-( const CartVect3D &a, const CartVect3D &b )
74 { return CartVect3D(a.x-b.x, a.y-b.y, a.z-b.z); }
75
operator %(const CartVect3D & a,const CartVect3D & b)76 static double operator%( const CartVect3D &a, const CartVect3D &b )
77 { return a.x*b.x + a.y*b.y + a.z*b.z; }
78
operator *(const CartVect3D & a,const CartVect3D & b)79 static CartVect3D operator*( const CartVect3D &a, const CartVect3D &b )
80 {
81 CartVect3D result;
82 result.x = a.y * b.z - a.z * b.y;
83 result.y = a.z * b.x - a.x * b.z;
84 result.z = a.x * b.y - a.y * b.x;
85 return result;
86 }
87
operator *=(const CartVect3D & o)88 CartVect3D& CartVect3D::operator*=( const CartVect3D& o )
89 { *this = *this * o; return *this; }
90
find_rotation(CartVect3D plane_normal,double matrix[3][3])91 static void find_rotation( CartVect3D plane_normal,
92 double matrix[3][3] )
93 {
94 // normalize
95 plane_normal /= plane_normal.len();
96 if (fabs(plane_normal.x) < 0.1)
97 plane_normal.x = 0.0;
98 if (fabs(plane_normal.y) < 0.1)
99 plane_normal.y = 0.0;
100 if (fabs(plane_normal.z) < 0.1)
101 plane_normal.z = 0.0;
102
103 // calculate vector to rotate about
104 const CartVect3D Z(0,0,1);
105 CartVect3D vector = plane_normal * Z;
106 const double len = vector.len();
107
108 // If vector is zero, no rotation
109 if (len < 1e-2) {
110 matrix[0][0] = matrix[1][1] = matrix[2][2] = 1.0;
111 matrix[0][1] = matrix[1][0] = 0.0;
112 matrix[0][2] = matrix[2][0] = 0.0;
113 matrix[1][2] = matrix[2][1] = 0.0;
114 return;
115 }
116 vector /= len;
117
118 const double cosine = plane_normal % Z;
119 const double sine = sqrt( 1 - cosine*cosine );
120
121 std::cerr << "Rotation: " << acos(cosine) << " [" << vector.x << ' ' << vector.y << ' ' << vector.z << ']' << std::endl;
122
123 const double x = vector.x;
124 const double y = vector.y;
125 const double z = vector.z;
126 const double c = cosine;
127 const double s = sine;
128 const double o = 1.0 - cosine;
129
130 matrix[0][0] = c + x*x*o;
131 matrix[0][1] = -z*s + x*y*o;
132 matrix[0][2] = y*s + x*z*o;
133
134 matrix[1][0] = z*s + x*z*o;
135 matrix[1][1] = c + y*y*o;
136 matrix[1][2] = -x*s + y*z*o;
137
138 matrix[2][0] = -y*s + x*z*o;
139 matrix[2][1] = x*s + y*z*o;
140 matrix[2][2] = c + z*z*o;
141 }
142
transform_point(CartVect3D & point,double matrix[3][3])143 static void transform_point( CartVect3D& point, double matrix[3][3] )
144 {
145 const double x = point.x;
146 const double y = point.y;
147 const double z = point.z;
148
149 point.x = x*matrix[0][0] + y*matrix[0][1] + z*matrix[0][2];
150 point.y = x*matrix[1][0] + y*matrix[1][1] + z*matrix[1][2];
151 point.z = x*matrix[2][0] + y*matrix[2][1] + z*matrix[2][2];
152 }
153
154
155 static void write_gnuplot( std::ostream& stream,
156 const std::vector<CartVect3D>& list );
157
158 static void write_svg( std::ostream& stream,
159 const std::vector<CartVect3D>& list );
160
161 static void write_eps( std::ostream& stream,
162 const std::vector<CartVect3D>& list,
163 int surface_id );
164
165 enum FileType {
166 POSTSCRIPT,
167 GNUPLOT,
168 SVG
169 };
170
171 using namespace moab;
172
main(int argc,char * argv[])173 int main(int argc, char* argv[])
174 {
175 Interface* moab = new Core();
176 ErrorCode result;
177 std::vector<CartVect3D>::iterator iter;
178 FileType type = GNUPLOT;
179
180 int idx = 1;
181 if (argc == 4)
182 {
183 if (!strcmp(argv[idx],"-p"))
184 type = POSTSCRIPT;
185 else if (!strcmp(argv[idx],"-g"))
186 type = GNUPLOT;
187 else if (!strcmp(argv[idx],"-s"))
188 type = SVG;
189 else
190 usage_error(argv[0]);
191 ++idx;
192 }
193
194
195
196 // scan CL args
197 int surface_id;
198 if (argc - idx != 2)
199 usage_error(argv[0]);
200 char* endptr;
201 surface_id = strtol( argv[idx], &endptr, 0 );
202 if (!endptr || *endptr)
203 usage_error(argv[0]);
204 ++idx;
205
206 // Load mesh
207 result = moab->load_mesh( argv[idx] );
208 if (MB_SUCCESS != result) {
209 if (MB_FILE_DOES_NOT_EXIST == result)
210 std::cerr << argv[idx] << " : open failed.\n";
211 else
212 std::cerr << argv[idx] << " : error reading file.\n";
213 return READ_ERROR;
214 }
215
216 // Get tag handles
217 EntityHandle surface;
218 const int dimension = 2; // surface
219 if (surface_id) {
220 Tag tags[2];
221 result = moab->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, tags[0] );
222 if (MB_SUCCESS != result) {
223 std::cerr << "No geometry tag.\n";
224 return OTHER_ERROR;
225 }
226 result = moab->tag_get_handle( GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, tags[1] );
227 if (MB_SUCCESS != result) {
228 std::cerr << "No ID tag.\n";
229 return OTHER_ERROR;
230 }
231
232 // Find entityset for surface.
233 const void* tag_values[] = { &dimension, &surface_id };
234 Range surfaces;
235 moab->get_entities_by_type_and_tag( 0, MBENTITYSET,
236 tags, tag_values,
237 2, surfaces );
238 if (surfaces.size() != 1) {
239 std::cerr << "Found " << surfaces.size()
240 << " surfaces with ID " << surface_id
241 << std::endl;
242 return SURFACE_NOT_FOUND;
243 }
244 surface = *surfaces.begin();
245 }
246 else {
247 surface = 0;
248 }
249
250 // Get surface mesh
251 Range elements;
252 result = moab->get_entities_by_dimension( surface, dimension, elements );
253 if (MB_SUCCESS != result) {
254 std::cerr << "Internal error\n";
255 return OTHER_ERROR;
256 }
257
258 // Calculate average corner normal in surface mesh
259 CartVect3D normal(0,0,0);
260 std::vector<EntityHandle> vertices;
261 std::vector<CartVect3D> coords;
262 for (Range::iterator i = elements.begin(); i != elements.end(); ++i)
263 {
264 vertices.clear();
265 result = moab->get_connectivity( &*i, 1, vertices, true );
266 if (MB_SUCCESS != result) {
267 std::cerr << "Internal error\n";
268 return OTHER_ERROR;
269 }
270 coords.clear();
271 coords.resize( vertices.size() );
272 result = moab->get_coords( &vertices[0], vertices.size(),
273 reinterpret_cast<double*>(&coords[0]) );
274 if (MB_SUCCESS != result) {
275 std::cerr << "Internal error\n";
276 return OTHER_ERROR;
277 }
278
279 for (size_t j = 0; j < coords.size(); ++j) {
280 CartVect3D v1 = coords[(j+1)%coords.size()] - coords[j];
281 CartVect3D v2 = coords[(j+1)%coords.size()] - coords[(j+2)%coords.size()];
282 normal += (v1 * v2);
283 }
284 }
285 normal /= normal.len();
286
287
288 // Get edges from elements
289 Range edge_range;
290 result = moab->get_adjacencies( elements, 1, true, edge_range, Interface::UNION );
291 if (MB_SUCCESS != result) {
292 std::cerr << "Internal error\n";
293 return OTHER_ERROR;
294 }
295
296 // Get vertex coordinates for each edge
297 std::vector<EntityHandle> edges( edge_range.size() );
298 std::copy( edge_range.begin(), edge_range.end(), edges.begin() );
299 vertices.clear();
300 result = moab->get_connectivity( &edges[0], edges.size(), vertices, true );
301 if (MB_SUCCESS != result) {
302 std::cerr << "Internal error\n";
303 return OTHER_ERROR;
304 }
305 coords.clear();
306 coords.resize( vertices.size() );
307 result = moab->get_coords( &vertices[0], vertices.size(),
308 reinterpret_cast<double*>(&coords[0]) );
309 if (MB_SUCCESS != result) {
310 std::cerr << "Internal error\n";
311 return OTHER_ERROR;
312 }
313
314 // Rotate points such that the projection into the view plane
315 // can be accomplished by discarding the 'z' coordinate of each
316 // point.
317
318 std::cerr << "Plane normal: [" << normal.x << ' ' << normal.y << ' ' << normal.z << ']' << std::endl;
319 double transform[3][3];
320 find_rotation( normal, transform );
321
322 for (iter = coords.begin(); iter != coords.end(); ++iter)
323 transform_point( *iter, transform );
324
325 // Write the file.
326
327 switch (type) {
328 case POSTSCRIPT:
329 write_eps( std::cout, coords, surface_id );
330 break;
331 case SVG:
332 write_svg( std::cout, coords );
333 break;
334 default:
335 write_gnuplot( std::cout, coords );
336 break;
337 }
338 return 0;
339 }
340
write_gnuplot(std::ostream & stream,const std::vector<CartVect3D> & coords)341 void write_gnuplot( std::ostream& stream,
342 const std::vector<CartVect3D>& coords )
343 {
344 std::vector<CartVect3D>::const_iterator iter;
345
346 stream << std::endl;
347 for (iter = coords.begin(); iter != coords.end(); ++iter) {
348 stream << iter->x << ' ' << iter->y << std::endl;
349 ++iter;
350 if (iter == coords.end()) {
351 stream << std::endl;
352 break;
353 }
354 stream << iter->x << ' ' << iter->y << std::endl;
355 stream << std::endl;
356 }
357 std::cerr << "Display with gnuplot command \"plot with lines\"\n";
358 }
359
box_max(CartVect3D & cur_max,const CartVect3D & pt)360 static void box_max( CartVect3D& cur_max, const CartVect3D& pt )
361 {
362 if (pt.x > cur_max.x)
363 cur_max.x = pt.x;
364 if (pt.y > cur_max.y)
365 cur_max.y = pt.y;
366 //if (pt.z > cur_max.z)
367 // cur_max.z = pt.z;
368 }
369
box_min(CartVect3D & cur_min,const CartVect3D & pt)370 static void box_min( CartVect3D& cur_min, const CartVect3D& pt )
371 {
372 if (pt.x < cur_min.x)
373 cur_min.x = pt.x;
374 if (pt.y < cur_min.y)
375 cur_min.y = pt.y;
376 //if (pt.z > cur_min.z)
377 // cur_min.z = pt.z;
378 }
379
380
write_eps(std::ostream & s,const std::vector<CartVect3D> & coords,int id)381 void write_eps( std::ostream& s, const std::vector<CartVect3D>& coords, int id )
382 {
383 // Coordinate range to use within EPS file
384 const int X_MAX = 540; // 540 pts / 72 pts/inch = 7.5 inches
385 const int Y_MAX = 720; // 720 pts / 72 pts/inch = 10 inches
386
387 std::vector<CartVect3D>::const_iterator iter;
388
389 // Get bounding box
390 const double D_MAX = std::numeric_limits<double>::max();
391 CartVect3D min( D_MAX, D_MAX, 0 );
392 CartVect3D max( -D_MAX, -D_MAX, 0 );
393 for (iter = coords.begin(); iter != coords.end(); ++iter)
394 {
395 box_max( max, *iter );
396 box_min( min, *iter );
397 }
398
399 // Calculate translation to page coordinates
400 CartVect3D offset = CartVect3D(0,0,0) - min;
401 CartVect3D scale = max - min;
402 scale.x = X_MAX / scale.x;
403 scale.y = Y_MAX / scale.y;
404 if (scale.x > scale.y) // keep proportions
405 scale.x = scale.y;
406 else
407 scale.y = scale.x;
408
409 //std::cerr << "Min: " << min.x << ' ' << min.y <<
410 // " Max: " << max.x << ' ' << max.y << std::endl
411 // << "Offset: " << offset.x << ' ' << offset.y <<
412 // " Scale: " << scale.x << ' ' << scale.y << std::endl;
413
414 // Write header stuff
415 s << "%!PS-Adobe-2.0 EPSF-2.0" << std::endl;
416 s << "%%Creator: MOAB surfplot" << std::endl;
417 s << "%%Title: Surface " << id << std::endl;
418 s << "%%DocumentData: Clean7Bit" << std::endl;
419 s << "%%Origin: 0 0" << std::endl;
420 int max_x = (int)((max.x + offset.x) * scale.x);
421 int max_y = (int)((max.y + offset.y) * scale.y);
422 s << "%%BoundingBox: 0 0 " << max_x << ' ' << max_y << std::endl;
423 s << "%%Pages: 1" << std::endl;
424
425 s << "%%BeginProlog" << std::endl;
426 s << "save" << std::endl;
427 s << "countdictstack" << std::endl;
428 s << "mark" << std::endl;
429 s << "newpath" << std::endl;
430 s << "/showpage {} def" << std::endl;
431 s << "/setpagedevice {pop} def" << std::endl;
432 s << "%%EndProlog" << std::endl;
433
434 s << "%%Page: 1 1" << std::endl;
435 s << "1 setlinewidth" << std::endl;
436 s << "0.0 setgray" << std::endl;
437
438 for (iter = coords.begin(); iter != coords.end(); ++iter)
439 {
440 double x1 = (iter->x + offset.x) * scale.x;
441 double y1 = (iter->y + offset.y) * scale.y;
442 if (++iter == coords.end())
443 break;
444 double x2 = (iter->x + offset.x) * scale.x;
445 double y2 = (iter->y + offset.y) * scale.y;
446
447 s << "newpath" << std::endl;
448 s << x1 << ' ' << y1 << " moveto" << std::endl;
449 s << x2 << ' ' << y2 << " lineto" << std::endl;
450 s << "stroke" << std::endl;
451 }
452
453 s << "%%Trailer" << std::endl;
454 s << "cleartomark" << std::endl;
455 s << "countdictstack" << std::endl;
456 s << "exch sub { end } repeat" << std::endl;
457 s << "restore" << std::endl;
458 s << "%%EOF" << std::endl;
459
460 }
461
462
write_svg(std::ostream & file,const std::vector<CartVect3D> & coords)463 void write_svg( std::ostream& file,
464 const std::vector<CartVect3D>& coords )
465 {
466
467 std::vector<CartVect3D>::const_iterator iter;
468
469 // Get bounding box
470 const double D_MAX = std::numeric_limits<double>::max();
471 CartVect3D min( D_MAX, D_MAX, 0 );
472 CartVect3D max( -D_MAX, -D_MAX, 0 );
473 for (iter = coords.begin(); iter != coords.end(); ++iter)
474 {
475 box_max( max, *iter );
476 box_min( min, *iter );
477 }
478 CartVect3D size = max - min;
479
480 // scale to 640 pixels on a side
481 double scale = 640.0 / (size.x > size.y ? size.x : size.y);
482 size *= scale;
483
484
485
486 file << "<?xml version=\"1.0\" standalone=\"no\"?>" << std::endl;
487 file << "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" "
488 << "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << std::endl;
489 file << std::endl;
490 file << "<svg width=\"" << (int)size.x
491 << "\" height=\"" << (int)size.y
492 << "\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" << std::endl;
493
494 int left = (int)(min.x * scale);
495 int top = (int)(min.y * scale);
496 iter = coords.begin();
497 while (iter != coords.end()) {
498 file << "<line "
499 << "x1=\"" << (int)(scale * iter->x) - left << "\" "
500 << "y1=\"" << (int)(scale * iter->y) - top << "\" ";
501 ++iter;
502 file << "x2=\"" << (int)(scale * iter->x) - left << "\" "
503 << "y2=\"" << (int)(scale * iter->y) - top << "\" "
504 << " style=\"stroke:rgb(99,99,99);stroke-width:2\""
505 << "/>" << std::endl;
506 ++iter;
507 }
508
509 // Write footer
510 file << "</svg>" << std::endl;
511 }
512
513
514
515
516
517
518