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