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 (const double *A, const double *B, const double *C, double Radius)
double	area_spherical_triangle (const double *A, const double *B, const double *C, double Radius)
double	area_spherical_polygon (const 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 (const double *ptA, const double *ptB, const double *ptC, double Radius)
double	area_spherical_polygon_lHuiller (const double *A, int N, double Radius, int *sign=NULL)
double	area_spherical_triangle_girard (const double *A, const double *B, const double *C, double Radius)
double	area_spherical_polygon_girard (const double *A, int N, double Radius)

Private Attributes
AreaMethod	m_eAreaMethod

Detailed Description

Definition at line 208 of file IntxUtils.hpp.

Member Enumeration Documentation

AreaMethod

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

IntxAreaUtils()

moab::IntxAreaUtils::IntxAreaUtils ( AreaMethod  p_eAreaMethod = lHuiller )

inline

Definition at line 218 of file IntxUtils.hpp.

: m_eAreaMethod( p_eAreaMethod ) {}

~IntxAreaUtils()

moab::IntxAreaUtils::~IntxAreaUtils ( )

inline

Definition at line 220 of file IntxUtils.hpp.

{}

Member Function Documentation

area_on_sphere()

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

Definition at line 1201 of file IntxUtils.cpp.

{
 // 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.
 if( elem_area <= 0 )
 {
 std::cout << "Area of element " << mb->id_from_handle( eh ) << " is = " << elem_area << "\n";
 mb->list_entity( eh );
 }
 assert( elem_area > 0 );
 
 // sum up the contribution
 total_area += elem_area;
 }
 
 // return total mesh area
 return total_area;
}

References area_spherical_element(), moab::Range::begin(), moab::Range::end(), ErrorCode, mb, MB_CHK_ERR_RET_VAL, MB_SUCCESS, moab::R, and moab::Range::size().

Referenced by main().

area_spherical_element()

double moab::IntxAreaUtils::area_spherical_element	(	Interface *	mb,
		EntityHandle	elem,
		double	R
	)

Definition at line 1244 of file IntxUtils.cpp.

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

References area_spherical_polygon(), ErrorCode, mb, MB_CHK_ERR_RET_VAL, and moab::R.

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

area_spherical_polygon()

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

Definition at line 911 of file IntxUtils.cpp.

{
 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 * Radius; //area_spherical_polygon_GQ normalizes
#endif
 case lHuiller:
 default:
 return area_spherical_polygon_lHuiller( A, N, Radius, sign );
 }
}

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

Referenced by area_spherical_element(), and main().

area_spherical_polygon_girard()

double moab::IntxAreaUtils::area_spherical_polygon_girard ( const double *  A, int  N, double  Radius  )

private

Definition at line 941 of file IntxUtils.cpp.

{
 // 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;
}

References moab::IntxUtils::oriented_spherical_angle().

Referenced by area_spherical_polygon().

area_spherical_polygon_lHuiller()

double moab::IntxAreaUtils::area_spherical_polygon_lHuiller ( const double *  A, int  N, double  Radius, int *  sign = NULL  )

private

Definition at line 958 of file IntxUtils.cpp.

{
 // 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;
}

References area_spherical_triangle_lHuiller().

Referenced by area_spherical_polygon(), and positive_orientation().

area_spherical_triangle()

double moab::IntxAreaUtils::area_spherical_triangle	(	const double *	A,
		const double *	B,
		const double *	C,
		double	Radius
	)

Definition at line 895 of file IntxUtils.cpp.

{
 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 );
#endif
 case lHuiller:
 default:
 return area_spherical_triangle_lHuiller( A, B, C, Radius );
 }
}

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

area_spherical_triangle_girard()

double moab::IntxAreaUtils::area_spherical_triangle_girard ( const double *  A, const double *  B, const double *  C, double  Radius  )

private

Definition at line 927 of file IntxUtils.cpp.

{
 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;
}

References spherical_angle().

Referenced by area_spherical_triangle().

area_spherical_triangle_lHuiller()

double moab::IntxAreaUtils::area_spherical_triangle_lHuiller ( const double *  ptA, const double *  ptB, const double *  ptC, double  Radius  )

private

Definition at line 1145 of file IntxUtils.cpp.

{
 
 // now, a is angle BOC, O is origin
 CartVect vA( ptA ), vB( ptB ), vC( ptC );
 double a = angle_robust( vB, vC );
 double b = angle_robust( vC, vA );
 double c = angle_robust( vA, vB );
 int sign = 1;
 // if( fabs( ( vA * vB ) % vC ) < 1e-17 ) sign = -1;
 if( ( vA * vB ) % vC < 0 ) sign = -1;
 double s = ( a + b + c ) / 2;
 double a1 = ( s - a ) / 2;
 double b1 = ( s - b ) / 2;
 double c1 = ( s - c ) / 2;
#ifdef MOAB_HAVE_TEMPESTREMAP
 if( fabs( a1 ) < 1.e-14 || fabs( b1 ) < 1.e-14 || fabs( c1 ) < 1.e-14 )
 {
 double area = area_spherical_triangle_GQ( ptA, ptB, ptC ) * sign;
#ifdef VERBOSE
 std::cout << " very obtuse angle, use TR to compute area " << " a1:" << a1 << " b1:" << b1 << " c1:" << c1
 << "\n";
 std::cout << " area with TR: " << area << "\n";
#endif
 return area;
 }
#endif
 double tmp = tan( s / 2 ) * tan( a1 ) * tan( b1 ) * tan( c1 );
 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;
}

References moab::angle_robust(), and moab::E.

Referenced by area_spherical_polygon_lHuiller(), area_spherical_triangle(), and positive_orientation().

positive_orientation()

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

Definition at line 1462 of file IntxUtils.cpp.

{
 Range cells2d;
 ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );MB_CHK_ERR( 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 );MB_CHK_ERR( rval );
 if( num_nodes < 3 ) return MB_FAILURE;
 
 double coords[9];
 rval = mb->get_coords( conn, 3, coords );MB_CHK_ERR( 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] );MB_CHK_ERR( rval );
 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 );MB_CHK_ERR( rval );
 }
 else
 {
 std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
 }
 }
 }
 return MB_SUCCESS;
}

References moab::IntxUtils::area2D(), area_spherical_polygon_lHuiller(), area_spherical_triangle_lHuiller(), moab::Range::begin(), moab::Range::end(), ErrorCode, mb, MB_CHK_ERR, MB_SUCCESS, and moab::R.

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

spherical_angle()

double moab::IntxAreaUtils::spherical_angle	(	const double *	A,
		const double *	B,
		const double *	C,
		double	Radius
	)

Definition at line 857 of file IntxUtils.cpp.

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

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

Referenced by area_spherical_triangle_girard().

Member Data Documentation

m_eAreaMethod

AreaMethod moab::IntxAreaUtils::m_eAreaMethod

private

Definition at line 268 of file IntxUtils.hpp.

Referenced by area_spherical_polygon(), and area_spherical_triangle().

---

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

• IntxUtils.hpp
• IntxUtils.cpp