moab
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   if (1 == argc) {
00024     std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
00025     return 0;
00026   }
00027   
00028   // declare variables
00029   Tag gtag, idtag;
00030   ErrorCode rval;
00031   const char *tag_geom = "GEOM_DIMENSION";
00032   const char *tag_gid = "GLOBAL_ID";
00033   Range sets;
00034   std::vector<EntityHandle> ents;
00035   
00036   // load a file
00037   mb = new Core();
00038   rval = mb->load_file(argv[1]);
00039   if (MB_SUCCESS != rval) return 1;
00040 
00041   // get the tag handle for the tags
00042   rval = mb->tag_get_handle(tag_geom, 1, MB_TYPE_INTEGER, gtag);
00043   if (MB_SUCCESS != rval) return 1;
00044   rval = mb->tag_get_handle(tag_gid, 1, MB_TYPE_INTEGER, idtag);
00045   if (MB_SUCCESS != rval) return 1;
00046 
00047   // get all the sets with GEOM_DIMESION tag
00048   rval = mb->get_entities_by_type_and_tag(0, MBENTITYSET, &gtag,
00049                       NULL, 1, sets);
00050   if (MB_SUCCESS != rval) return 1;
00051  
00052   // iterate over each set, getting entities
00053   Range::iterator set_it;
00054     
00055   // loop thru all the geometric entity sets
00056   for (set_it = sets.begin(); set_it != sets.end(); set_it++)  {
00057 
00058     EntityHandle this_set = *set_it;
00059     
00060     // get the id for this set
00061     int set_id;
00062     rval = mb->tag_get_data(gtag, &this_set, 1, &set_id);
00063     if (MB_SUCCESS != rval) return 1;
00064 
00065     // check if it is a surface entities (GEOM_DIMENSION 2) then compute area
00066     if (set_id ==2){
00067       
00068       // area of a surface
00069       double total_area=0.0;
00070 
00071       //get the global id of this surface
00072       int gid = 0;
00073       rval = mb->tag_get_data(idtag, &this_set, 1, &gid);
00074       if (MB_SUCCESS != rval) return 1;
00075 
00076       // get all entities with dimension 2 in ents
00077       rval = mb->get_entities_by_dimension(this_set, 2, ents);
00078       if (MB_SUCCESS != rval) return 1;
00079  
00080       // compute the area
00081       total_area = compute_area(ents);
00082 
00083       ents.clear();
00084       
00085       std::cout << "Total area of meshes in surface " << gid << " =  " << total_area << std::endl;
00086     }
00087   }
00088 }
00089 
00090 // This routine takes all the element entities in a face as input and computes the surface area
00091 // iterating over each element
00092 double compute_area(std::vector<EntityHandle> & entities){
00093  
00094   ErrorCode rval= MB_SUCCESS;
00095   double area = 0.0;
00096  
00097   // loop thro' all the elements
00098   for (int i=0;i<int(entities.size());i++){
00099     std::vector<EntityHandle> conn;
00100     EntityHandle handle = entities[i];
00101 
00102     // get the connectivity of this element
00103     rval = mb->get_connectivity(&handle, 1, conn);
00104     if (MB_SUCCESS != rval) return -1.0;
00105 
00106     // break polygon into triangles and sum the area - Limitation: Convex polygon
00107     for (int j = 2; j<=int(conn.size()); ++j){
00108 
00109       EntityHandle vertices[3]={conn[0], conn[j-1],  conn[j-2]};
00110       CartVect coords[3];
00111       
00112       // get 3 coordinates forming the triangle
00113       rval = mb->get_coords(vertices, 3, coords[0].array());
00114       if (MB_SUCCESS != rval) return -1.0;
00115       
00116       CartVect edge0 = coords[1] - coords [0];
00117       CartVect edge1 = coords [2] - coords[0];
00118       
00119       // using MBCarVect overloaded operators and computing triangle area
00120       area+=(edge0*edge1).length()/2.0;
00121     }
00122   }
00123   // clear the entities, else old entities remain
00124   entities.clear();
00125   return area;
00126 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines