MOAB: Mesh Oriented datABase  (version 5.2.1)
moab::IntxAreaUtils Class Reference

#include <IntxUtils.hpp>

Public Types

enum  AreaMethod { lHuiller = 0, Girard = 1, GaussQuadrature = 2 }

Public Member Functions

 IntxAreaUtils (AreaMethod p_eAreaMethod=lHuiller)
 ~IntxAreaUtils ()
double spherical_angle (double *A, double *B, double *C, double Radius)
double area_spherical_triangle (double *A, double *B, double *C, double Radius)
double area_spherical_polygon (double *A, int N, double Radius, int *sign=NULL)
double area_spherical_element (Interface *mb, EntityHandle elem, double R)
double area_on_sphere (Interface *mb, EntityHandle set, double R)
ErrorCode positive_orientation (Interface *mb, EntityHandle set, double R)

Private Member Functions

double area_spherical_triangle_lHuiller (double *ptA, double *ptB, double *ptC, double Radius)
double area_spherical_polygon_lHuiller (double *A, int N, double Radius, int *sign=NULL)
double area_spherical_triangle_girard (double *A, double *B, double *C, double Radius)
double area_spherical_polygon_girard (double *A, int N, double Radius)

Private Attributes

AreaMethod m_eAreaMethod

Detailed Description

Definition at line 170 of file IntxUtils.hpp.


Member Enumeration Documentation

Enumerator:
lHuiller 
Girard 
GaussQuadrature 

Definition at line 173 of file IntxUtils.hpp.

    {
        lHuiller        = 0,
        Girard          = 1,
        GaussQuadrature = 2
    };

Constructor & Destructor Documentation

Definition at line 180 of file IntxUtils.hpp.

: m_eAreaMethod( p_eAreaMethod ) {}

Definition at line 182 of file IntxUtils.hpp.

{}

Member Function Documentation

double moab::IntxAreaUtils::area_on_sphere ( Interface *  mb,
EntityHandle  set,
double  R 
)

Definition at line 982 of file IntxUtils.cpp.

References area_spherical_element(), moab::Range::begin(), moab::Range::end(), ie, MB_CHK_ERR_RET_VAL, MB_SUCCESS, and moab::Range::size().

Referenced by main(), test_intx_in_parallel_elem_based(), and test_intx_mpas().

{
    // Get all entities of dimension 2
    Range inputRange;
    ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );

    // Filter by elements that are owned by current process
    std::vector< int > ownerinfo( inputRange.size(), -1 );
    Tag intxOwnerTag;
    rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
    if( MB_SUCCESS == rval )
    {
        rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
    }

    // compare total area with 4*M_PI * R^2
    int ie            = 0;
    double total_area = 0.;
    for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
    {

        // All zero/positive owner data represents ghosted elems
        if( ownerinfo[ie++] >= 0 ) continue;

        EntityHandle eh        = *eit;
        const double elem_area = this->area_spherical_element( mb, eh, R );

        // check whether the area of the spherical element is positive.
        assert( elem_area > 0 );

        // sum up the contribution
        total_area += elem_area;
    }

    // return total mesh area
    return total_area;
}
double moab::IntxAreaUtils::area_spherical_element ( Interface *  mb,
EntityHandle  elem,
double  R 
)

Definition at line 1020 of file IntxUtils.cpp.

References area_spherical_polygon(), and MB_CHK_ERR_RET_VAL.

Referenced by area_on_sphere(), and moab::Intx2MeshOnSphere::update_tracer_data().

{
    // get the nodes, then the coordinates
    const EntityHandle* verts;
    int nsides;
    ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );

    // account for possible padded polygons
    while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
        nsides--;

    // get coordinates
    std::vector< double > coords( 3 * nsides );
    rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );

    // compute and return the area of the polygonal element
    return area_spherical_polygon( &coords[0], nsides, R );
}
double moab::IntxAreaUtils::area_spherical_polygon ( double *  A,
int  N,
double  Radius,
int *  sign = NULL 
)

Definition at line 813 of file IntxUtils.cpp.

References area_spherical_polygon_girard(), area_spherical_polygon_lHuiller(), GaussQuadrature, Girard, lHuiller, and m_eAreaMethod.

Referenced by add_field_value(), area_spherical_element(), and main().

{
    switch( m_eAreaMethod )
    {
        case Girard:
            return area_spherical_polygon_girard( A, N, Radius );
#ifdef MOAB_HAVE_TEMPESTREMAP
        case GaussQuadrature:
            return area_spherical_polygon_GQ( A, N, Radius );
#endif
        case lHuiller:
        default:
            return area_spherical_polygon_lHuiller( A, N, Radius, sign );
    }
}
double moab::IntxAreaUtils::area_spherical_polygon_girard ( double *  A,
int  N,
double  Radius 
) [private]

Definition at line 843 of file IntxUtils.cpp.

References N, and moab::IntxUtils::oriented_spherical_angle().

Referenced by area_spherical_polygon().

{
    // this should work for non-convex polygons too
    // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
    //
    if( N <= 2 ) return 0.;
    double sum_angles = 0.;
    for( int i = 0; i < N; i++ )
    {
        int i1 = ( i + 1 ) % N;
        int i2 = ( i + 2 ) % N;
        sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
    }
    double correction = sum_angles - ( N - 2 ) * M_PI;
    return Radius * Radius * correction;
}
double moab::IntxAreaUtils::area_spherical_polygon_lHuiller ( double *  A,
int  N,
double  Radius,
int *  sign = NULL 
) [private]

Definition at line 860 of file IntxUtils.cpp.

References area_spherical_triangle_lHuiller().

Referenced by area_spherical_polygon().

{
    // This should work for non-convex polygons too
    // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
    // We also assume that the orientation is positive;
    // If negative orientation, the area will be negative
    if( N <= 2 ) return 0.;

    int lsign   = 1;  // assume positive orientain
    double area = 0.;
    for( int i = 1; i < N - 1; i++ )
    {
        int i1              = i + 1;
        double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
        if( areaTriangle < 0 )
            lsign = -1;  // signal that we have at least one triangle with negative orientation ;
                         // possible nonconvex polygon
        area += areaTriangle;
    }
    if( sign ) *sign = lsign;

    return area;
}
double moab::IntxAreaUtils::area_spherical_triangle ( double *  A,
double *  B,
double *  C,
double  Radius 
)

Definition at line 797 of file IntxUtils.cpp.

References area_spherical_triangle_girard(), area_spherical_triangle_lHuiller(), GaussQuadrature, Girard, lHuiller, and m_eAreaMethod.

{
    switch( m_eAreaMethod )
    {
        case Girard:
            return area_spherical_triangle_girard( A, B, C, Radius );
#ifdef MOAB_HAVE_TEMPESTREMAP
        case GaussQuadrature:
            return area_spherical_triangle_GQ( A, B, C, Radius );
#endif
        case lHuiller:
        default:
            return area_spherical_triangle_lHuiller( A, B, C, Radius );
    }
}
double moab::IntxAreaUtils::area_spherical_triangle_girard ( double *  A,
double *  B,
double *  C,
double  Radius 
) [private]

Definition at line 829 of file IntxUtils.cpp.

References b, and spherical_angle().

Referenced by area_spherical_triangle().

{
    double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
                        spherical_angle( C, A, B, Radius ) - M_PI;
    double area = Radius * Radius * correction;
    // now, is it negative or positive? is it pointing toward the center or outward?
    CartVect a( A ), b( B ), c( C );
    CartVect abc = ( b - a ) * ( c - a );
    if( abc % a > 0 )  // dot product positive, means ABC points out
        return area;
    else
        return -area;
}
double moab::IntxAreaUtils::area_spherical_triangle_lHuiller ( double *  ptA,
double *  ptB,
double *  ptC,
double  Radius 
) [private]

Definition at line 944 of file IntxUtils.cpp.

References moab::angle(), b, moab::E, Radius, and MBMesquite::sign.

Referenced by area_spherical_polygon_lHuiller(), and area_spherical_triangle().

{

    // now, a is angle BOC, O is origin
    CartVect vA( ptA ), vB( ptB ), vC( ptC );
    double a = angle( vB, vC );
    double b = angle( vC, vA );
    double c = angle( vA, vB );
    int sign = 1;
    if( ( vA * vB ) % vC < 0 ) sign = -1;
    double s   = ( a + b + c ) / 2;
    double tmp = tan( s / 2 ) * tan( ( s - a ) / 2 ) * tan( ( s - b ) / 2 ) * tan( ( s - c ) / 2 );
    if( tmp < 0. ) tmp = 0.;

    double E = 4 * atan( sqrt( tmp ) );
    if( E != E ) std::cout << " NaN at spherical triangle area \n";

    double area = sign * E * Radius * Radius;

#ifdef CHECKNEGATIVEAREA
    if( area < 0 )
    {
        std::cout << "negative area: " << area << "\n";
        std::cout << std::setprecision( 15 );
        std::cout << "vA: " << vA << "\n";
        std::cout << "vB: " << vB << "\n";
        std::cout << "vC: " << vC << "\n";
        std::cout << "sign: " << sign << "\n";
        std::cout << " a: " << a << "\n";
        std::cout << " b: " << b << "\n";
        std::cout << " c: " << c << "\n";
    }
#endif

    return area;
}
ErrorCode moab::IntxAreaUtils::positive_orientation ( Interface *  mb,
EntityHandle  set,
double  R 
)

Definition at line 1236 of file IntxUtils.cpp.

References moab::IntxUtils::area2D(), moab::Range::begin(), conn, moab::Range::end(), MB_FAILURE, and MB_SUCCESS.

Referenced by moab::TempestRemapper::ComputeOverlapMesh(), and main().

{
    Range cells2d;
    ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );
    if( MB_SUCCESS != rval ) return rval;
    for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
    {
        EntityHandle cell        = *qit;
        const EntityHandle* conn = NULL;
        int num_nodes            = 0;
        rval                     = mb->get_connectivity( cell, conn, num_nodes );
        if( MB_SUCCESS != rval ) return rval;
        if( num_nodes < 3 ) return MB_FAILURE;

        double coords[9];
        rval = mb->get_coords( conn, 3, coords );
        if( MB_SUCCESS != rval ) return rval;

        double area;
        if( R > 0 )
            area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
        else
            area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
        if( area < 0 )
        {
            // compute all area, do not revert if total area is positive
            std::vector< double > coords2( 3 * num_nodes );
            // get coordinates
            rval = mb->get_coords( conn, num_nodes, &coords2[0] );
            if( MB_SUCCESS != rval ) return MB_FAILURE;
            double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
            if( totArea < 0 )
            {
                std::vector< EntityHandle > newconn( num_nodes );
                for( int i = 0; i < num_nodes; i++ )
                {
                    newconn[num_nodes - 1 - i] = conn[i];
                }
                rval = mb->set_connectivity( cell, &newconn[0], num_nodes );
                if( MB_SUCCESS != rval ) return rval;
            }
            else
            {
                std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
            }
        }
    }
    return MB_SUCCESS;
}
double moab::IntxAreaUtils::spherical_angle ( double *  A,
double *  B,
double *  C,
double  Radius 
)

Definition at line 761 of file IntxUtils.cpp.

References moab::angle(), b, moab::CartVect::length_squared(), and Radius.

Referenced by area_spherical_triangle_girard().

{
    // the angle by definition is between the planes OAB and OBC
    CartVect a( A );
    CartVect b( B );
    CartVect c( C );
    double err1 = a.length_squared() - Radius * Radius;
    if( fabs( err1 ) > 0.0001 )
    { std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n"; }
    CartVect normalOAB = a * b;
    CartVect normalOCB = c * b;
    return angle( normalOAB, normalOCB );
}

Member Data Documentation

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines