#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.

218 : m_eAreaMethod( p_eAreaMethod ) {}

◆ ~IntxAreaUtils()

moab::IntxAreaUtils::~IntxAreaUtils ( )

inline

Definition at line 220 of file IntxUtils.hpp.

220 {}

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:

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ AreaMethod

Constructor & Destructor Documentation

◆ IntxAreaUtils()

◆ ~IntxAreaUtils()

Member Function Documentation

◆ area_on_sphere()

◆ area_spherical_element()

◆ area_spherical_polygon()

◆ area_spherical_polygon_girard()

◆ area_spherical_polygon_lHuiller()

◆ area_spherical_triangle()

◆ area_spherical_triangle_girard()

◆ area_spherical_triangle_lHuiller()

◆ positive_orientation()

◆ spherical_angle()

Member Data Documentation

◆ m_eAreaMethod