moab
|
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, >ag, 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 }