Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
SurfArea.cpp
Go to the documentation of this file.
1 /*
2  This example takes a .cub mesh file as input and prints out the total area of
3  meshes in each surface of the model. It works for tri and quad elements.
4  It makes use of CUBIT's reserved tag - GEOM_DIMENSION and GLOBAL_ID tag.
5  Both GLOBAL_ID & GEOM_DIMENSION tag are associated with all the geometric
6  entities in a .cub file. Note: The program would give incorrect result for a
7  non-convex element, since it breaks polygons into triangles for computing the area
8 */
9 
10 #include "moab/Core.hpp"
11 #include "moab/Range.hpp"
12 #include "moab/CartVect.hpp"
13 #include <iostream>
14 
15 using namespace moab;
16 
17 double compute_area( std::vector< EntityHandle >& );
18 
19 // instantiate
21 
22 int main( int argc, char** argv )
23 {
24  if( 1 == argc )
25  {
26  std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
27  return 0;
28  }
29 
30  // declare variables
31  Tag gtag, idtag;
32  ErrorCode rval;
33  const char* tag_geom = "GEOM_DIMENSION";
34  const char* tag_gid = "GLOBAL_ID";
35  Range sets;
36  std::vector< EntityHandle > ents;
37 
38  // load a file
39  mb = new Core();
40  rval = mb->load_file( argv[1] );
41  if( MB_SUCCESS != rval ) return 1;
42 
43  // get the tag handle for the tags
44  rval = mb->tag_get_handle( tag_geom, 1, MB_TYPE_INTEGER, gtag );
45  if( MB_SUCCESS != rval ) return 1;
46  rval = mb->tag_get_handle( tag_gid, 1, MB_TYPE_INTEGER, idtag );
47  if( MB_SUCCESS != rval ) return 1;
48 
49  // get all the sets with GEOM_DIMESION tag
50  rval = mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &gtag, NULL, 1, sets );
51  if( MB_SUCCESS != rval ) return 1;
52 
53  // iterate over each set, getting entities
54  Range::iterator set_it;
55 
56  // loop thru all the geometric entity sets
57  for( set_it = sets.begin(); set_it != sets.end(); ++set_it )
58  {
59 
60  EntityHandle this_set = *set_it;
61 
62  // get the id for this set
63  int set_id;
64  rval = mb->tag_get_data( gtag, &this_set, 1, &set_id );
65  if( MB_SUCCESS != rval ) return 1;
66 
67  // check if it is a surface entities (GEOM_DIMENSION 2) then compute area
68  if( set_id == 2 )
69  {
70 
71  // area of a surface
72  double total_area = 0.0;
73 
74  // get the global id of this surface
75  int gid = 0;
76  rval = mb->tag_get_data( idtag, &this_set, 1, &gid );
77  if( MB_SUCCESS != rval ) return 1;
78 
79  // get all entities with dimension 2 in ents
80  rval = mb->get_entities_by_dimension( this_set, 2, ents );
81  if( MB_SUCCESS != rval ) return 1;
82 
83  // compute the area
84  total_area = compute_area( ents );
85 
86  ents.clear();
87 
88  std::cout << "Total area of meshes in surface " << gid << " = " << total_area << std::endl;
89  }
90  }
91 }
92 
93 // This routine takes all the element entities in a face as input and computes the surface area
94 // iterating over each element
95 double compute_area( std::vector< EntityHandle >& entities )
96 {
97 
98  ErrorCode rval = MB_SUCCESS;
99  double area = 0.0;
100 
101  // loop thro' all the elements
102  for( int i = 0; i < int( entities.size() ); i++ )
103  {
104  std::vector< EntityHandle > conn;
105  EntityHandle handle = entities[i];
106 
107  // get the connectivity of this element
108  rval = mb->get_connectivity( &handle, 1, conn );
109  if( MB_SUCCESS != rval ) return -1.0;
110 
111  // break polygon into triangles and sum the area - Limitation: Convex polygon
112  for( int j = 2; j <= int( conn.size() ); ++j )
113  {
114 
115  EntityHandle vertices[3] = { conn[0], conn[j - 1], conn[j - 2] };
116  CartVect coords[3];
117 
118  // get 3 coordinates forming the triangle
119  rval = mb->get_coords( vertices, 3, coords[0].array() );
120  if( MB_SUCCESS != rval ) return -1.0;
121 
122  CartVect edge0 = coords[1] - coords[0];
123  CartVect edge1 = coords[2] - coords[0];
124 
125  // using MBCarVect overloaded operators and computing triangle area
126  area += ( edge0 * edge1 ).length() / 2.0;
127  }
128  }
129  // clear the entities, else old entities remain
130  entities.clear();
131  return area;
132 }