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
   20 Interface* mb;
   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 }