 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#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 |
Definition at line 208 of file IntxUtils.hpp.
Definition at line 211 of file IntxUtils.hpp.
{
lHuiller = 0,
Girard = 1,
GaussQuadrature = 2
};
| 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.
{}
| 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 );
}
AreaMethod moab::IntxAreaUtils::m_eAreaMethod [private] |
Definition at line 268 of file IntxUtils.hpp.
Referenced by area_spherical_polygon(), and area_spherical_triangle().
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.