#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, int rank=0)
double	area_spherical_polygon (double A, int N, double Radius, int sign=NULL, int rank=0)
double	area_spherical_element (Interface *mb, EntityHandle elem, double R, int rank=0)
double	area_on_sphere (Interface *mb, EntityHandle set, double R, int rank=0)
ErrorCode	positive_orientation (Interface *mb, EntityHandle set, double R, int rank=0)
Private Member Functions
double	area_spherical_triangle_lHuiller (double ptA, double ptB, double *ptC, double Radius, int rank=0)
double	area_spherical_polygon_lHuiller (double A, int N, double Radius, int sign=NULL, int rank=0)
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 208 of file IntxUtils.hpp.

Member Enumeration Documentation

enum moab::IntxAreaUtils::AreaMethod

Enumerator:

lHuiller
Girard
GaussQuadrature

Definition at line 211 of file IntxUtils.hpp.

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

Constructor & Destructor Documentation

moab::IntxAreaUtils::IntxAreaUtils ( AreaMethod p_eAreaMethod = lHuiller ) [inline]

Definition at line 218 of file IntxUtils.hpp.

: m_eAreaMethod( p_eAreaMethod ) {}

moab::IntxAreaUtils::~IntxAreaUtils ( ) [inline]

Definition at line 220 of file IntxUtils.hpp.

{}

Member Function Documentation

double moab::IntxAreaUtils::area_on_sphere	(	Interface *	mb,
		EntityHandle	set,
		double	R,
		int	rank = `0`
	)

Definition at line 1018 of file IntxUtils.cpp.

References area_spherical_element(), moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_entities_by_dimension(), MB_CHK_ERR_RET_VAL, MB_SUCCESS, moab::Range::size(), moab::Interface::tag_get_data(), and moab::Interface::tag_get_handle().

Referenced by main(), and test_intx_in_parallel_elem_based().

{
    // 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, rank );

        // 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,
		int	rank = `0`
	)

Definition at line 1056 of file IntxUtils.cpp.

References area_spherical_polygon(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::list_entity(), 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
    int lsign   = 1;
    double area = area_spherical_polygon( &coords[0], nsides, R, &lsign, rank );
    {
        if( lsign < 0 )
        {
            std::cout << " IntxAreaUtils::area_spherical_element : lsign :" << lsign << " area:" << area << " rank:"
                    << rank << "\n";
            mb->list_entity( elem );
        }
    }
    return area;
}

double moab::IntxAreaUtils::area_spherical_polygon	(	double *	A,
		int	N,
		double	Radius,
		int *	sign = `NULL`,
		int	rank = `0`
	)

Definition at line 849 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, rank );
    }
}

double moab::IntxAreaUtils::area_spherical_polygon_girard	(	double *	A,
		int	N,
		double	Radius
	)		`[private]`

Definition at line 879 of file IntxUtils.cpp.

References 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`,
		int	rank = `0`
	)		`[private]`

Definition at line 896 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 orientation
    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, rank );
        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,
		int	rank = `0`
	)

Definition at line 833 of file IntxUtils.cpp.

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

Referenced by moab::Intx2MeshOnSphere::findNodes().

{
    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, rank );
    }
}

double moab::IntxAreaUtils::area_spherical_triangle_girard	(	double *	A,
		double *	B,
		double *	C,
		double	Radius
	)		`[private]`

Definition at line 865 of file IntxUtils.cpp.

References 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,
		int	rank = `0`
	)		`[private]`

Definition at line 980 of file IntxUtils.cpp.

References moab::angle(), moab::E, and Radius.

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 << " on task :" << rank << "\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,
		int	rank = `0`
	)

Definition at line 1282 of file IntxUtils.cpp.

References moab::IntxUtils::area2D(), moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_dimension(), moab::Interface::globalId_tag(), MB_SUCCESS, moab::Interface::set_connectivity(), and moab::Interface::tag_get_data().

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;
    Tag gidTag = mb->globalId_tag();
    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, rank );
        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;
            int lsign;
            double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R, &lsign, rank );
            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;
                int gid;
                rval = mb->tag_get_data( gidTag, &cell, 1, &gid );
                if( MB_SUCCESS != rval ) return rval;
                std::cout << " revert cell with global id:" << gid << " with num_nodes " << num_nodes << " on rank "
                          << rank << "\n";
            }
            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 795 of file IntxUtils.cpp.

References moab::angle(), 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

AreaMethod moab::IntxAreaUtils::m_eAreaMethod [private]

Definition at line 268 of file IntxUtils.hpp.

Referenced by area_spherical_polygon(), and area_spherical_triangle().

List of all members.

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation