SurfArea.cpp

Go to the documentation of this file.

/* 
   This example takes a .cub mesh file as input and prints out the total area of 
   meshes in each surface of the model. It works for tri and quad elements.
   It makes use of CUBIT's reserved tag - GEOM_DIMENSION and GLOBAL_ID tag.
   Both GLOBAL_ID & GEOM_DIMENSION tag are associated with all the geometric 
   entities in a .cub file. Note: The program would give incorrect result for a 
   non-convex element, since it breaks polygons into triangles for computing the area
*/

#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/CartVect.hpp"
#include <iostream>

using namespace moab;

double compute_area(std::vector<EntityHandle>&);

// instantiate
Interface *mb;

int main(int argc, char **argv) {
  if (1 == argc) {
    std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
    return 0;
  }
  
  // declare variables
  Tag gtag, idtag;
  ErrorCode rval;
  const char *tag_geom = "GEOM_DIMENSION";
  const char *tag_gid = "GLOBAL_ID";
  Range sets;
  std::vector<EntityHandle> ents;
  
  // load a file
  mb = new Core();
  rval = mb->load_file(argv[1]);
  if (MB_SUCCESS != rval) return 1;

  // get the tag handle for the tags
  rval = mb->tag_get_handle(tag_geom, 1, MB_TYPE_INTEGER, gtag);
  if (MB_SUCCESS != rval) return 1;
  rval = mb->tag_get_handle(tag_gid, 1, MB_TYPE_INTEGER, idtag);
  if (MB_SUCCESS != rval) return 1;

  // get all the sets with GEOM_DIMESION tag
  rval = mb->get_entities_by_type_and_tag(0, MBENTITYSET, &gtag,
                      NULL, 1, sets);
  if (MB_SUCCESS != rval) return 1;
 
  // iterate over each set, getting entities
  Range::iterator set_it;
    
  // loop thru all the geometric entity sets
  for (set_it = sets.begin(); set_it != sets.end(); set_it++)  {

    EntityHandle this_set = *set_it;
    
    // get the id for this set
    int set_id;
    rval = mb->tag_get_data(gtag, &this_set, 1, &set_id);
    if (MB_SUCCESS != rval) return 1;

    // check if it is a surface entities (GEOM_DIMENSION 2) then compute area
    if (set_id ==2){
      
      // area of a surface
      double total_area=0.0;

      //get the global id of this surface
      int gid = 0;
      rval = mb->tag_get_data(idtag, &this_set, 1, &gid);
      if (MB_SUCCESS != rval) return 1;

      // get all entities with dimension 2 in ents
      rval = mb->get_entities_by_dimension(this_set, 2, ents);
      if (MB_SUCCESS != rval) return 1;
 
      // compute the area
      total_area = compute_area(ents);

      ents.clear();
      
      std::cout << "Total area of meshes in surface " << gid << " =  " << total_area << std::endl;
    }
  }
}

// This routine takes all the element entities in a face as input and computes the surface area
// iterating over each element
double compute_area(std::vector<EntityHandle> & entities){
 
  ErrorCode rval= MB_SUCCESS;
  double area = 0.0;
 
  // loop thro' all the elements
  for (int i=0;i<int(entities.size());i++){
    std::vector<EntityHandle> conn;
    EntityHandle handle = entities[i];

    // get the connectivity of this element
    rval = mb->get_connectivity(&handle, 1, conn);
    if (MB_SUCCESS != rval) return -1.0;

    // break polygon into triangles and sum the area - Limitation: Convex polygon
    for (int j = 2; j<=int(conn.size()); ++j){

      EntityHandle vertices[3]={conn[0], conn[j-1],  conn[j-2]};
      CartVect coords[3];
      
      // get 3 coordinates forming the triangle
      rval = mb->get_coords(vertices, 3, coords[0].array());
      if (MB_SUCCESS != rval) return -1.0;
      
      CartVect edge0 = coords[1] - coords [0];
      CartVect edge1 = coords [2] - coords[0];
      
      // using MBCarVect overloaded operators and computing triangle area
      area+=(edge0*edge1).length()/2.0;
    }
  }
  // clear the entities, else old entities remain
  entities.clear();
  return area;
}