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 (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

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

1202 {
1203  // Get all entities of dimension 2
1204  Range inputRange;
1205  ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1206 
1207  // Filter by elements that are owned by current process
1208  std::vector< int > ownerinfo( inputRange.size(), -1 );
1209  Tag intxOwnerTag;
1210  rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
1211  if( MB_SUCCESS == rval )
1212  {
1213  rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1214  }
1215 
1216  // compare total area with 4*M_PI * R^2
1217  int ie = 0;
1218  double total_area = 0.;
1219  for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
1220  {
1221 
1222  // All zero/positive owner data represents ghosted elems
1223  if( ownerinfo[ie++] >= 0 ) continue;
1224 
1225  EntityHandle eh = *eit;
1226  const double elem_area = this->area_spherical_element( mb, eh, R );
1227 
1228  // check whether the area of the spherical element is positive.
1229  if( elem_area <= 0 )
1230  {
1231  std::cout << "Area of element " << mb->id_from_handle( eh ) << " is = " << elem_area << "\n";
1232  mb->list_entity( eh );
1233  }
1234  assert( elem_area > 0 );
1235 
1236  // sum up the contribution
1237  total_area += elem_area;
1238  }
1239 
1240  // return total mesh area
1241  return total_area;
1242 }

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.

1245 {
1246  // get the nodes, then the coordinates
1247  const EntityHandle* verts;
1248  int nsides;
1249  ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1250 
1251  // account for possible padded polygons
1252  while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1253  nsides--;
1254 
1255  // get coordinates
1256  std::vector< double > coords( 3 * nsides );
1257  rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
1258 
1259  // compute and return the area of the polygonal element
1260  return area_spherical_polygon( &coords[0], nsides, R );
1261 }

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.

912 {
913  switch( m_eAreaMethod )
914  {
915  case Girard:
916  return area_spherical_polygon_girard( A, N, Radius );
917 #ifdef MOAB_HAVE_TEMPESTREMAP
918  case GaussQuadrature:
919  return area_spherical_polygon_GQ( A, N ) * Radius * Radius; //area_spherical_polygon_GQ normalizes
920 #endif
921  case lHuiller:
922  default:
923  return area_spherical_polygon_lHuiller( A, N, Radius, sign );
924  }
925 }

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.

942 {
943  // this should work for non-convex polygons too
944  // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
945  //
946  if( N <= 2 ) return 0.;
947  double sum_angles = 0.;
948  for( int i = 0; i < N; i++ )
949  {
950  int i1 = ( i + 1 ) % N;
951  int i2 = ( i + 2 ) % N;
952  sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
953  }
954  double correction = sum_angles - ( N - 2 ) * M_PI;
955  return Radius * Radius * correction;
956 }

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.

959 {
960  // This should work for non-convex polygons too
961  // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
962  // We also assume that the orientation is positive;
963  // If negative orientation, the area will be negative
964  if( N <= 2 ) return 0.;
965 
966  int lsign = 1; // assume positive orientain
967  double area = 0.;
968  for( int i = 1; i < N - 1; i++ )
969  {
970  int i1 = i + 1;
971  double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
972  if( areaTriangle < 0 )
973  lsign = -1; // signal that we have at least one triangle with negative orientation ;
974  // possible nonconvex polygon
975  area += areaTriangle;
976  }
977  if( sign ) *sign = lsign;
978 
979  return area;
980 }

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.

896 {
897  switch( m_eAreaMethod )
898  {
899  case Girard:
900  return area_spherical_triangle_girard( A, B, C, Radius );
901 #ifdef MOAB_HAVE_TEMPESTREMAP
902  case GaussQuadrature:
903  return area_spherical_triangle_GQ( A, B, C );
904 #endif
905  case lHuiller:
906  default:
907  return area_spherical_triangle_lHuiller( A, B, C, Radius );
908  }
909 }

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.

928 {
929  double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
930  spherical_angle( C, A, B, Radius ) - M_PI;
931  double area = Radius * Radius * correction;
932  // now, is it negative or positive? is it pointing toward the center or outward?
933  CartVect a( A ), b( B ), c( C );
934  CartVect abc = ( b - a ) * ( c - a );
935  if( abc % a > 0 ) // dot product positive, means ABC points out
936  return area;
937  else
938  return -area;
939 }

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.

1149 {
1150 
1151  // now, a is angle BOC, O is origin
1152  CartVect vA( ptA ), vB( ptB ), vC( ptC );
1153  double a = angle_robust( vB, vC );
1154  double b = angle_robust( vC, vA );
1155  double c = angle_robust( vA, vB );
1156  int sign = 1;
1157  // if( fabs( ( vA * vB ) % vC ) < 1e-17 ) sign = -1;
1158  if( ( vA * vB ) % vC < 0 ) sign = -1;
1159  double s = ( a + b + c ) / 2;
1160  double a1 = ( s - a ) / 2;
1161  double b1 = ( s - b ) / 2;
1162  double c1 = ( s - c ) / 2;
1163 #ifdef MOAB_HAVE_TEMPESTREMAP
1164  if( fabs( a1 ) < 1.e-14 || fabs( b1 ) < 1.e-14 || fabs( c1 ) < 1.e-14 )
1165  {
1166  double area = area_spherical_triangle_GQ( ptA, ptB, ptC ) * sign;
1167 #ifdef VERBOSE
1168  std::cout << " very obtuse angle, use TR to compute area " << " a1:" << a1 << " b1:" << b1 << " c1:" << c1
1169  << "\n";
1170  std::cout << " area with TR: " << area << "\n";
1171 #endif
1172  return area;
1173  }
1174 #endif
1175  double tmp = tan( s / 2 ) * tan( a1 ) * tan( b1 ) * tan( c1 );
1176  if( tmp < 0. ) tmp = 0.;
1177 
1178  double E = 4 * atan( sqrt( tmp ) );
1179  if( E != E ) std::cout << " NaN at spherical triangle area \n";
1180 
1181  double area = sign * E * Radius * Radius;
1182 
1183 #ifdef CHECKNEGATIVEAREA
1184  if( area < 0 )
1185  {
1186  std::cout << "negative area: " << area << "\n";
1187  std::cout << std::setprecision( 15 );
1188  std::cout << "vA: " << vA << "\n";
1189  std::cout << "vB: " << vB << "\n";
1190  std::cout << "vC: " << vC << "\n";
1191  std::cout << "sign: " << sign << "\n";
1192  std::cout << " a: " << a << "\n";
1193  std::cout << " b: " << b << "\n";
1194  std::cout << " c: " << c << "\n";
1195  }
1196 #endif
1197 
1198  return area;
1199 }

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.

1463 {
1464  Range cells2d;
1465  ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );MB_CHK_ERR( rval );
1466  for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
1467  {
1468  EntityHandle cell = *qit;
1469  const EntityHandle* conn = NULL;
1470  int num_nodes = 0;
1471  rval = mb->get_connectivity( cell, conn, num_nodes );MB_CHK_ERR( rval );
1472  if( num_nodes < 3 ) return MB_FAILURE;
1473 
1474  double coords[9];
1475  rval = mb->get_coords( conn, 3, coords );MB_CHK_ERR( rval );
1476 
1477  double area;
1478  if( R > 0 )
1479  area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
1480  else
1481  area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
1482  if( area < 0 )
1483  {
1484  // compute all area, do not revert if total area is positive
1485  std::vector< double > coords2( 3 * num_nodes );
1486  // get coordinates
1487  rval = mb->get_coords( conn, num_nodes, &coords2[0] );MB_CHK_ERR( rval );
1488  double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
1489  if( totArea < 0 )
1490  {
1491  std::vector< EntityHandle > newconn( num_nodes );
1492  for( int i = 0; i < num_nodes; i++ )
1493  {
1494  newconn[num_nodes - 1 - i] = conn[i];
1495  }
1496  rval = mb->set_connectivity( cell, &newconn[0], num_nodes );MB_CHK_ERR( rval );
1497  }
1498  else
1499  {
1500  std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
1501  }
1502  }
1503  }
1504  return MB_SUCCESS;
1505 }

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.

858 {
859  // the angle by definition is between the planes OAB and OBC
860  CartVect a( A );
861  CartVect b( B );
862  CartVect c( C );
863  double err1 = a.length_squared() - Radius * Radius;
864  if( fabs( err1 ) > 0.0001 )
865  {
866  std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n";
867  }
868  CartVect normalOAB = a * b;
869  CartVect normalOCB = c * b;
870  return angle( normalOAB, normalOCB );
871 }

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: