MOAB: Mesh Oriented datABase  (version 5.4.1)
SurfArea.cpp
Go to the documentation of this file.
00001 /*
00002    This example takes a .cub mesh file as input and prints out the total area of
00003    meshes in each surface of the model. It works for tri and quad elements.
00004    It makes use of CUBIT's reserved tag - GEOM_DIMENSION and GLOBAL_ID tag.
00005    Both GLOBAL_ID & GEOM_DIMENSION tag are associated with all the geometric
00006    entities in a .cub file. Note: The program would give incorrect result for a
00007    non-convex element, since it breaks polygons into triangles for computing the area
00008 */
00009 
00010 #include "moab/Core.hpp"
00011 #include "moab/Range.hpp"
00012 #include "moab/CartVect.hpp"
00013 #include <iostream>
00014 
00015 using namespace moab;
00016 
00017 double compute_area( std::vector< EntityHandle >& );
00018 
00019 // instantiate
00020 Interface* mb;
00021 
00022 int main( int argc, char** argv )
00023 {
00024     if( 1 == argc )
00025     {
00026         std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
00027         return 0;
00028     }
00029 
00030     // declare variables
00031     Tag gtag, idtag;
00032     ErrorCode rval;
00033     const char* tag_geom = "GEOM_DIMENSION";
00034     const char* tag_gid  = "GLOBAL_ID";
00035     Range sets;
00036     std::vector< EntityHandle > ents;
00037 
00038     // load a file
00039     mb   = new Core();
00040     rval = mb->load_file( argv[1] );
00041     if( MB_SUCCESS != rval ) return 1;
00042 
00043     // get the tag handle for the tags
00044     rval = mb->tag_get_handle( tag_geom, 1, MB_TYPE_INTEGER, gtag );
00045     if( MB_SUCCESS != rval ) return 1;
00046     rval = mb->tag_get_handle( tag_gid, 1, MB_TYPE_INTEGER, idtag );
00047     if( MB_SUCCESS != rval ) return 1;
00048 
00049     // get all the sets with GEOM_DIMESION tag
00050     rval = mb->get_entities_by_type_and_tag( 0, MBENTITYSET, &gtag, NULL, 1, sets );
00051     if( MB_SUCCESS != rval ) return 1;
00052 
00053     // iterate over each set, getting entities
00054     Range::iterator set_it;
00055 
00056     // loop thru all the geometric entity sets
00057     for( set_it = sets.begin(); set_it != sets.end(); ++set_it )
00058     {
00059 
00060         EntityHandle this_set = *set_it;
00061 
00062         // get the id for this set
00063         int set_id;
00064         rval = mb->tag_get_data( gtag, &this_set, 1, &set_id );
00065         if( MB_SUCCESS != rval ) return 1;
00066 
00067         // check if it is a surface entities (GEOM_DIMENSION 2) then compute area
00068         if( set_id == 2 )
00069         {
00070 
00071             // area of a surface
00072             double total_area = 0.0;
00073 
00074             // get the global id of this surface
00075             int gid = 0;
00076             rval    = mb->tag_get_data( idtag, &this_set, 1, &gid );
00077             if( MB_SUCCESS != rval ) return 1;
00078 
00079             // get all entities with dimension 2 in ents
00080             rval = mb->get_entities_by_dimension( this_set, 2, ents );
00081             if( MB_SUCCESS != rval ) return 1;
00082 
00083             // compute the area
00084             total_area = compute_area( ents );
00085 
00086             ents.clear();
00087 
00088             std::cout << "Total area of meshes in surface " << gid << " =  " << total_area << std::endl;
00089         }
00090     }
00091 }
00092 
00093 // This routine takes all the element entities in a face as input and computes the surface area
00094 // iterating over each element
00095 double compute_area( std::vector< EntityHandle >& entities )
00096 {
00097 
00098     ErrorCode rval = MB_SUCCESS;
00099     double area    = 0.0;
00100 
00101     // loop thro' all the elements
00102     for( int i = 0; i < int( entities.size() ); i++ )
00103     {
00104         std::vector< EntityHandle > conn;
00105         EntityHandle handle = entities[i];
00106 
00107         // get the connectivity of this element
00108         rval = mb->get_connectivity( &handle, 1, conn );
00109         if( MB_SUCCESS != rval ) return -1.0;
00110 
00111         // break polygon into triangles and sum the area - Limitation: Convex polygon
00112         for( int j = 2; j <= int( conn.size() ); ++j )
00113         {
00114 
00115             EntityHandle vertices[3] = { conn[0], conn[j - 1], conn[j - 2] };
00116             CartVect coords[3];
00117 
00118             // get 3 coordinates forming the triangle
00119             rval = mb->get_coords( vertices, 3, coords[0].array() );
00120             if( MB_SUCCESS != rval ) return -1.0;
00121 
00122             CartVect edge0 = coords[1] - coords[0];
00123             CartVect edge1 = coords[2] - coords[0];
00124 
00125             // using MBCarVect overloaded operators and computing triangle area
00126             area += ( edge0 * edge1 ).length() / 2.0;
00127         }
00128     }
00129     // clear the entities, else old entities remain
00130     entities.clear();
00131     return area;
00132 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines