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 

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] << " " << 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, >ag, 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;
}