Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 208 of file IntxUtils.hpp.

Member Enumeration Documentation

◆ AreaMethod

Enumerator
lHuiller 
Girard 
GaussQuadrature 

Definition at line 211 of file IntxUtils.hpp.

212  {
213  lHuiller = 0,
214  Girard = 1,
215  GaussQuadrature = 2
216  };

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 1018 of file IntxUtils.cpp.

1019 {
1020  // Get all entities of dimension 2
1021  Range inputRange;
1022  ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1023 
1024  // Filter by elements that are owned by current process
1025  std::vector< int > ownerinfo( inputRange.size(), -1 );
1026  Tag intxOwnerTag;
1027  rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
1028  if( MB_SUCCESS == rval )
1029  {
1030  rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1031  }
1032 
1033  // compare total area with 4*M_PI * R^2
1034  int ie = 0;
1035  double total_area = 0.;
1036  for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
1037  {
1038 
1039  // All zero/positive owner data represents ghosted elems
1040  if( ownerinfo[ie++] >= 0 ) continue;
1041 
1042  EntityHandle eh = *eit;
1043  const double elem_area = this->area_spherical_element( mb, eh, R );
1044 
1045  // check whether the area of the spherical element is positive.
1046  assert( elem_area > 0 );
1047 
1048  // sum up the contribution
1049  total_area += elem_area;
1050  }
1051 
1052  // return total mesh area
1053  return total_area;
1054 }

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 1056 of file IntxUtils.cpp.

1057 {
1058  // get the nodes, then the coordinates
1059  const EntityHandle* verts;
1060  int nsides;
1061  ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1062 
1063  // account for possible padded polygons
1064  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1065  nsides--;
1066 
1067  // get coordinates
1068  std::vector< double > coords( 3 * nsides );
1069  rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1070 
1071  // compute and return the area of the polygonal element
1072  return area_spherical_polygon( &coords[0], nsides, R );
1073 }

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 ( double *  A,
int  N,
double  Radius,
int *  sign = NULL 
)

Definition at line 849 of file IntxUtils.cpp.

850 {
851  switch( m_eAreaMethod )
852  {
853  case Girard:
854  return area_spherical_polygon_girard( A, N, Radius );
855 #ifdef MOAB_HAVE_TEMPESTREMAP
856  case GaussQuadrature:
857  return area_spherical_polygon_GQ( A, N, Radius );
858 #endif
859  case lHuiller:
860  default:
861  return area_spherical_polygon_lHuiller( A, N, Radius, sign );
862  }
863 }

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 ( double *  A,
int  N,
double  Radius 
)
private

Definition at line 879 of file IntxUtils.cpp.

880 {
881  // this should work for non-convex polygons too
882  // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
883  //
884  if( N <= 2 ) return 0.;
885  double sum_angles = 0.;
886  for( int i = 0; i < N; i++ )
887  {
888  int i1 = ( i + 1 ) % N;
889  int i2 = ( i + 2 ) % N;
890  sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
891  }
892  double correction = sum_angles - ( N - 2 ) * M_PI;
893  return Radius * Radius * correction;
894 }

References moab::IntxUtils::oriented_spherical_angle().

Referenced by area_spherical_polygon().

◆ area_spherical_polygon_lHuiller()

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

Definition at line 896 of file IntxUtils.cpp.

897 {
898  // This should work for non-convex polygons too
899  // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
900  // We also assume that the orientation is positive;
901  // If negative orientation, the area will be negative
902  if( N <= 2 ) return 0.;
903 
904  int lsign = 1; // assume positive orientain
905  double area = 0.;
906  for( int i = 1; i < N - 1; i++ )
907  {
908  int i1 = i + 1;
909  double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
910  if( areaTriangle < 0 )
911  lsign = -1; // signal that we have at least one triangle with negative orientation ;
912  // possible nonconvex polygon
913  area += areaTriangle;
914  }
915  if( sign ) *sign = lsign;
916 
917  return area;
918 }

References area_spherical_triangle_lHuiller().

Referenced by area_spherical_polygon(), and positive_orientation().

◆ area_spherical_triangle()

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

Definition at line 833 of file IntxUtils.cpp.

834 {
835  switch( m_eAreaMethod )
836  {
837  case Girard:
838  return area_spherical_triangle_girard( A, B, C, Radius );
839 #ifdef MOAB_HAVE_TEMPESTREMAP
840  case GaussQuadrature:
841  return area_spherical_triangle_GQ( A, B, C, Radius );
842 #endif
843  case lHuiller:
844  default:
845  return area_spherical_triangle_lHuiller( A, B, C, Radius );
846  }
847 }

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 ( double *  A,
double *  B,
double *  C,
double  Radius 
)
private

Definition at line 865 of file IntxUtils.cpp.

866 {
867  double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
868  spherical_angle( C, A, B, Radius ) - M_PI;
869  double area = Radius * Radius * correction;
870  // now, is it negative or positive? is it pointing toward the center or outward?
871  CartVect a( A ), b( B ), c( C );
872  CartVect abc = ( b - a ) * ( c - a );
873  if( abc % a > 0 ) // dot product positive, means ABC points out
874  return area;
875  else
876  return -area;
877 }

References spherical_angle().

Referenced by area_spherical_triangle().

◆ area_spherical_triangle_lHuiller()

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

Definition at line 980 of file IntxUtils.cpp.

981 {
982 
983  // now, a is angle BOC, O is origin
984  CartVect vA( ptA ), vB( ptB ), vC( ptC );
985  double a = angle( vB, vC );
986  double b = angle( vC, vA );
987  double c = angle( vA, vB );
988  int sign = 1;
989  if( ( vA * vB ) % vC < 0 ) sign = -1;
990  double s = ( a + b + c ) / 2;
991  double tmp = tan( s / 2 ) * tan( ( s - a ) / 2 ) * tan( ( s - b ) / 2 ) * tan( ( s - c ) / 2 );
992  if( tmp < 0. ) tmp = 0.;
993 
994  double E = 4 * atan( sqrt( tmp ) );
995  if( E != E ) std::cout << " NaN at spherical triangle area \n";
996 
997  double area = sign * E * Radius * Radius;
998 
999 #ifdef CHECKNEGATIVEAREA
1000  if( area < 0 )
1001  {
1002  std::cout << "negative area: " << area << "\n";
1003  std::cout << std::setprecision( 15 );
1004  std::cout << "vA: " << vA << "\n";
1005  std::cout << "vB: " << vB << "\n";
1006  std::cout << "vC: " << vC << "\n";
1007  std::cout << "sign: " << sign << "\n";
1008  std::cout << " a: " << a << "\n";
1009  std::cout << " b: " << b << "\n";
1010  std::cout << " c: " << c << "\n";
1011  }
1012 #endif
1013 
1014  return area;
1015 }

References moab::angle(), 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 1272 of file IntxUtils.cpp.

1273 {
1274  Range cells2d;
1275  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );
1276  if( MB_SUCCESS != rval ) return rval;
1277  for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
1278  {
1279  EntityHandle cell = *qit;
1280  const EntityHandle* conn = NULL;
1281  int num_nodes = 0;
1282  rval = mb->get_connectivity( cell, conn, num_nodes );
1283  if( MB_SUCCESS != rval ) return rval;
1284  if( num_nodes < 3 ) return MB_FAILURE;
1285 
1286  double coords[9];
1287  rval = mb->get_coords( conn, 3, coords );
1288  if( MB_SUCCESS != rval ) return rval;
1289 
1290  double area;
1291  if( R > 0 )
1292  area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
1293  else
1294  area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
1295  if( area < 0 )
1296  {
1297  // compute all area, do not revert if total area is positive
1298  std::vector< double > coords2( 3 * num_nodes );
1299  // get coordinates
1300  rval = mb->get_coords( conn, num_nodes, &coords2[0] );
1301  if( MB_SUCCESS != rval ) return MB_FAILURE;
1302  double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
1303  if( totArea < 0 )
1304  {
1305  std::vector< EntityHandle > newconn( num_nodes );
1306  for( int i = 0; i < num_nodes; i++ )
1307  {
1308  newconn[num_nodes - 1 - i] = conn[i];
1309  }
1310  rval = mb->set_connectivity( cell, &newconn[0], num_nodes );
1311  if( MB_SUCCESS != rval ) return rval;
1312  }
1313  else
1314  {
1315  std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
1316  }
1317  }
1318  }
1319  return MB_SUCCESS;
1320 }

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

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

◆ spherical_angle()

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

Definition at line 795 of file IntxUtils.cpp.

796 {
797  // the angle by definition is between the planes OAB and OBC
798  CartVect a( A );
799  CartVect b( B );
800  CartVect c( C );
801  double err1 = a.length_squared() - Radius * Radius;
802  if( fabs( err1 ) > 0.0001 )
803  {
804  std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n";
805  }
806  CartVect normalOAB = a * b;
807  CartVect normalOCB = c * b;
808  return angle( normalOAB, normalOCB );
809 }

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: