moab::IntxUtils Class Reference

#include <IntxUtils.hpp>

Classes
struct	SphereCoords

Static Public Member Functions
static double	dist2 (double *a, double *b)
static double	area2D (double *a, double *b, double *c)
static int	borderPointsOfXinY2 (double *X, int nX, double *Y, int nY, double *P, int *side, double epsilon_area)
static int	SortAndRemoveDoubles2 (double *P, int &nP, double epsilon)
static ErrorCode	EdgeIntersections2 (double *blue, int nsBlue, double *red, int nsRed, int *markb, int *markr, double *points, int &nPoints)
static ErrorCode	EdgeIntxRllCs (double *blue, CartVect *bluec, int *blueEdgeType, int nsBlue, double *red, CartVect *redc, int nsRed, int *markb, int *markr, int plane, double Radius, double *points, int &nPoints)
static void	decide_gnomonic_plane (const CartVect &pos, int &oPlane)
static ErrorCode	gnomonic_projection (const CartVect &pos, double R, int plane, double &c1, double &c2)
static ErrorCode	reverse_gnomonic_projection (const double &c1, const double &c2, double R, int plane, CartVect &pos)
static void	gnomonic_unroll (double &c1, double &c2, double R, int plane)
static ErrorCode	global_gnomonic_projection (Interface *mb, EntityHandle inSet, double R, bool centers_only, EntityHandle &outSet)
static void	transform_coordinates (double *avg_position, int projection_type)
static SphereCoords	cart_to_spherical (CartVect &)
static CartVect	spherical_to_cart (SphereCoords &)
static ErrorCode	ScaleToRadius (Interface *mb, EntityHandle set, double R)
static double	distance_on_great_circle (CartVect &p1, CartVect &p2)
static ErrorCode	enforce_convexity (Interface *mb, EntityHandle set, int rank=0)
	Enforces convexity for a given set of polygons. More...
static double	oriented_spherical_angle (const double *A, const double *B, const double *C)
static ErrorCode	fix_degenerate_quads (Interface *mb, EntityHandle set)
static double	distance_on_sphere (double la1, double te1, double la2, double te2)
static ErrorCode	intersect_great_circle_arcs (double *A, double *B, double *C, double *D, double R, double *E)
static ErrorCode	intersect_great_circle_arc_with_clat_arc (double *A, double *B, double *C, double *D, double R, double *E, int &np)
static int	borderPointsOfCSinRLL (CartVect *redc, double *red2dc, int nsRed, CartVect *bluec, int nsBlue, int *blueEdgeType, double *P, int *side, double epsil)
static ErrorCode	deep_copy_set_with_quads (Interface *mb, EntityHandle source_set, EntityHandle dest_set)
static ErrorCode	remove_duplicate_vertices (Interface *mb, EntityHandle file_set, double merge_tol, std::vector< Tag > &tagList)
static ErrorCode	remove_padded_vertices (Interface *mb, EntityHandle file_set, std::vector< Tag > &tagList)

Detailed Description

Definition at line 16 of file IntxUtils.hpp.

Member Function Documentation

area2D()

static double moab::IntxUtils::area2D ( double *  a, double *  b, double *  c  )

inlinestatic

Definition at line 26 of file IntxUtils.hpp.

 {
 // (b-a)x(c-a) / 2
 return ( ( b[0] - a[0] ) * ( c[1] - a[1] ) - ( b[1] - a[1] ) * ( c[0] - a[0] ) ) / 2;
 }

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), moab::IntxAreaUtils::positive_orientation(), moab::Intx2MeshInPlane::setup_tgt_cell(), moab::Intx2MeshOnSphere::setup_tgt_cell(), and moab::IntxRllCssphere::setup_tgt_cell().

borderPointsOfCSinRLL()

int moab::IntxUtils::borderPointsOfCSinRLL ( CartVect *  redc, double *  red2dc, int  nsRed, CartVect *  bluec, int  nsBlue, int *  blueEdgeType, double *  P, int *  side, double  epsil  )

static

Definition at line 1918 of file IntxUtils.cpp.

{
 int extraPoints = 0;
 // first decide the blue z coordinates
 CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
 for( int i = 0; i < nsBlue; i++ )
 {
 if( blueEdgeType[i] == 0 )
 {
 int iP1 = ( i + 1 ) % nsBlue;
 if( bluec[i][2] > bluec[iP1][2] )
 {
 A = bluec[i];
 B = bluec[iP1];
 C = bluec[( i + 2 ) % nsBlue];
 D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
 break;
 }
 }
 }
 if( nsBlue == 3 && B[2] < 0 )
 {
 // select D to be C
 D = C;
 C = B; // B is the south pole then
 }
 // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
 // AB is const longitude, BC and DA constant latitude
 // check now each of the red points if they are inside this rectangle
 for( int i = 0; i < nsRed; i++ )
 {
 CartVect& X = redc[i];
 if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
 // now decide if it is between the planes OAB and OCD
 if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
 {
 side[i] = 1; //
 // it means point X is in the rectangle that we want , on the sphere
 // pass the coords 2d
 P[extraPoints * 2] = red2dc[2 * i];
 P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
 extraPoints++;
 }
 }
 return extraPoints;
}

Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

borderPointsOfXinY2()

int moab::IntxUtils::borderPointsOfXinY2 ( double *  X, int  nX, double *  Y, int  nY, double *  P, int *  side, double  epsilon_area  )

static

Computes the border points of X in Y2.

Parameters

X	The array of points representing X.
nX	The number of points in X.
Y	The array of points representing Y.
nY	The number of points in Y.
P	The array to store the border points of X in Y2.
side	The array to store the side information for each point in X.
epsilon_area	The epsilon value for area comparison.

Returns

The number of extra points found.

Definition at line 67 of file IntxUtils.cpp.

{
 // 2 triangles, 3 corners, is the corner of X in Y?
 // Y must have a positive area
 /*
 */
 int extraPoint = 0;
 for( int i = 0; i < nX; i++ )
 {
 // compute double the area of all nY triangles formed by a side of Y and a corner of X; if
 // one is negative, stop (negative means it is outside; X and Y are all oriented such that
 // they are positive oriented;
 // if one area is negative, it means it is outside the convex region, for sure)
 double* A = X + 2 * i;
 
 int inside = 1;
 for( int j = 0; j < nY; j++ )
 {
 double* B = Y + 2 * j;
 int j1 = ( j + 1 ) % nY;
 double* C = Y + 2 * j1; // no copy of data
 
 double area2 = ( B[0] - A[0] ) * ( C[1] - A[1] ) - ( C[0] - A[0] ) * ( B[1] - A[1] );
 if( area2 < -epsilon_area )
 {
 inside = 0;
 break;
 }
 }
 if( inside )
 {
 side[i] = 1; // so vertex i of X is inside the convex region formed by Y
 // so side has nX dimension (first array)
 P[extraPoint * 2] = A[0];
 P[extraPoint * 2 + 1] = A[1];
 extraPoint++;
 }
 }
 return extraPoint;
}

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), and moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

cart_to_spherical()

IntxUtils::SphereCoords moab::IntxUtils::cart_to_spherical ( CartVect &  cart3d )

static

Definition at line 788 of file IntxUtils.cpp.

{
 SphereCoords res;
 res.R = cart3d.length();
 if( res.R < 0 )
 {
 res.lon = res.lat = 0.;
 return res;
 }
 res.lat = asin( cart3d[2] / res.R );
 res.lon = atan2( cart3d[1], cart3d[0] );
 if( res.lon < 0 ) res.lon += 2 * M_PI; // M_PI is defined in math.h? it seems to be true, although
 // there are some defines it depends on :(
 // #if defined __USE_BSD || defined __USE_XOPEN ???
 
 return res;
}

References moab::IntxUtils::SphereCoords::lat, moab::CartVect::length(), moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.

Referenced by distance_on_great_circle().

decide_gnomonic_plane()

void moab::IntxUtils::decide_gnomonic_plane ( const CartVect &  pos, int &  oPlane  )

static

Definition at line 411 of file IntxUtils.cpp.

{
 // decide plane, based on max x, y, z
 if( fabs( pos[0] ) < fabs( pos[1] ) )
 {
 if( fabs( pos[2] ) < fabs( pos[1] ) )
 {
 // pos[1] is biggest
 if( pos[1] > 0 )
 plane = 2;
 else
 plane = 4;
 }
 else
 {
 // pos[2] is biggest
 if( pos[2] < 0 )
 plane = 5;
 else
 plane = 6;
 }
 }
 else
 {
 if( fabs( pos[2] ) < fabs( pos[0] ) )
 {
 // pos[0] is the greatest
 if( pos[0] > 0 )
 plane = 1;
 else
 plane = 3;
 }
 else
 {
 // pos[2] is biggest
 if( pos[2] < 0 )
 plane = 5;
 else
 plane = 6;
 }
 }
 return;
}

Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), global_gnomonic_projection(), moab::Intx2MeshOnSphere::setup_tgt_cell(), moab::IntxRllCssphere::setup_tgt_cell(), transform_coordinates(), and moab::BoundBox::update_box_spherical_elem().

deep_copy_set_with_quads()

ErrorCode moab::IntxUtils::deep_copy_set_with_quads ( Interface *  mb, EntityHandle  source_set, EntityHandle  dest_set  )

static

Definition at line 1973 of file IntxUtils.cpp.

{
 ReadUtilIface* read_iface;
 ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
 // create the handle tag for the corresponding element / vertex
 
 EntityHandle dum = 0;
 Tag corrTag = 0; // it will be created here
 rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
 
 // give the same global id to new verts and cells created in the lagr(departure) mesh
 Tag gid = mb->globalId_tag();
 
 Range quads;
 rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
 
 Range connecVerts;
 rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
 
 std::map< EntityHandle, EntityHandle > newNodes;
 
 std::vector< double* > coords;
 EntityHandle start_vert, start_elem, *connect;
 int num_verts = connecVerts.size();
 rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
 if( MB_SUCCESS != rval ) return rval;
 // fill it up
 int i = 0;
 for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
 {
 EntityHandle oldV = *vit;
 CartVect posi;
 rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
 
 int global_id;
 rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
 EntityHandle new_vert = start_vert + i;
 // Cppcheck warning (false positive): variable coords is assigned a value that is never used
 coords[0][i] = posi[0];
 coords[1][i] = posi[1];
 coords[2][i] = posi[2];
 
 newNodes[oldV] = new_vert;
 // set also the correspondent tag :)
 rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
 
 // also the other side
 // need to check if we really need this; the new vertex will never need the old vertex
 // we have the global id which is the same
 rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
 // set the global id on the corresponding vertex the same as the initial vertex
 rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
 }
 // now create new quads in order (in a sequence)
 
 rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
 if( MB_SUCCESS != rval ) return rval;
 int ie = 0;
 for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
 {
 EntityHandle q = *it;
 int nnodes;
 const EntityHandle* conn;
 rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
 int global_id;
 rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
 
 for( int ii = 0; ii < nnodes; ii++ )
 {
 EntityHandle v1 = conn[ii];
 connect[4 * ie + ii] = newNodes[v1];
 }
 EntityHandle newElement = start_elem + ie;
 
 // set the corresponding tag; not sure we need this one, from old to new
 rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
 rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
 
 // set the global id
 rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
 
 rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
 }
 
 rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
 
 return MB_SUCCESS;
}

References moab::Range::begin(), CORRTAGNAME, moab::dum, moab::Range::end(), ErrorCode, moab::ReadUtilIface::get_element_connect(), moab::ReadUtilIface::get_node_coords(), mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_HANDLE, MBQUAD, moab::Range::size(), and moab::ReadUtilIface::update_adjacencies().

dist2()

static double moab::IntxUtils::dist2 ( double *  a, double *  b  )

inlinestatic

Definition at line 20 of file IntxUtils.hpp.

 {
 double abx = b[0] - a[0], aby = b[1] - a[1];
 return sqrt( abx * abx + aby * aby );
 }

Referenced by moab::Intx2MeshInPlane::findNodes(), moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and SortAndRemoveDoubles2().

distance_on_great_circle()

double moab::IntxUtils::distance_on_great_circle ( CartVect &  p1, CartVect &  p2  )

static

Definition at line 1263 of file IntxUtils.cpp.

{
 SphereCoords sph1 = cart_to_spherical( p1 );
 SphereCoords sph2 = cart_to_spherical( p2 );
 // radius should be the same
 return sph1.R *
 acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
}

References cart_to_spherical(), moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.

distance_on_sphere()

double moab::IntxUtils::distance_on_sphere ( double  la1, double  te1, double  la2, double  te2  )

static

Definition at line 1509 of file IntxUtils.cpp.

{
 return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
}

EdgeIntersections2()

ErrorCode moab::IntxUtils::EdgeIntersections2 ( double *  blue, int  nsBlue, double *  red, int  nsRed, int *  markb, int *  markr, double *  points, int &  nPoints  )

static

Computes the edge intersections of two elements.

Parameters

blue	The array of points representing the blue element.
nsBlue	The number of points in the blue element.
red	The array of points representing the red element.
nsRed	The number of points in the red element.
markb	The array to mark the intersecting edges of the blue element.
markr	The array to mark the intersecting edges of the red element.
points	The array to store the intersection points.
nPoints	The number of intersection points found.

Returns

The error code.

Definition at line 222 of file IntxUtils.cpp.

{
 /* EDGEINTERSECTIONS computes edge intersections of two elements
 [P,n]=EdgeIntersections(X,Y) computes for the two given elements * red
 and blue ( stored column wise )
 (point coordinates are stored column-wise, in counter clock
 order) the points P where their edges intersect. In addition,
 in n the indices of which neighbors of red are also intersecting
 with blue are given.
 */
 
 // points is an array with enough slots (24 * 2 doubles)
 nPoints = 0;
 for( int i = 0; i < MAXEDGES; i++ )
 {
 markb[i] = markr[i] = 0;
 }
 
 for( int i = 0; i < nsBlue; i++ )
 {
 for( int j = 0; j < nsRed; j++ )
 {
 double b[2];
 double a[2][2]; // 2*2
 int iPlus1 = ( i + 1 ) % nsBlue;
 int jPlus1 = ( j + 1 ) % nsRed;
 for( int k = 0; k < 2; k++ )
 {
 b[k] = red[2 * j + k] - blue[2 * i + k];
 // row k of a: a(k, 0), a(k, 1)
 a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
 a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
 }
 double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
 if( fabs( delta ) > 1.e-14 ) // this is close to machine epsilon
 {
 // not parallel
 double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
 double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
 if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
 {
 // the intersection is good
 for( int k = 0; k < 2; k++ )
 {
 points[2 * nPoints + k] = blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
 }
 markb[i] = 1; // so neighbor number i of blue will be considered too.
 markr[j] = 1; // this will be used in advancing red around blue quad
 nPoints++;
 }
 }
 // the case delta ~ 0. will be considered by the interior points logic
 }
 }
 return MB_SUCCESS;
}

References MAXEDGES, and MB_SUCCESS.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), and moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().

EdgeIntxRllCs()

ErrorCode moab::IntxUtils::EdgeIntxRllCs ( double *  blue, CartVect *  bluec, int *  blueEdgeType, int  nsBlue, double *  red, CartVect *  redc, int  nsRed, int *  markb, int *  markr, int  plane, double  R, double *  points, int &  nPoints  )

static

Computes the edge intersections between a RLL and CS quad.

Note: Special function.

Parameters

blue	The array of points representing the blue element.
bluec	The array of Cartesian coordinates for the blue element.
blueEdgeType	The array of edge types for the blue element.
nsBlue	The number of points in the blue element.
red	The array of points representing the red element.
redc	The array of Cartesian coordinates for the red element.
nsRed	The number of points in the red element.
markb	The array to mark the intersecting edges of the blue element.
markr	The array to mark the intersecting edges of the red element.
plane	The plane of intersection.
R	The radius of the sphere.
points	The array to store the intersection points.
nPoints	The number of intersection points found.

Returns

The error code.

Definition at line 306 of file IntxUtils.cpp.

{
 // if blue edge type is 1, intersect in 3d then project to 2d by gnomonic projection
 // everything else the same (except if there are 2 points resulting, which is rare)
 for( int i = 0; i < 4; i++ )
 { // always at most 4 , so maybe don't bother
 markb[i] = markr[i] = 0;
 }
 
 for( int i = 0; i < nsBlue; i++ )
 {
 int iPlus1 = ( i + 1 ) % nsBlue;
 if( blueEdgeType[i] == 0 ) // old style, just 2d
 {
 for( int j = 0; j < nsRed; j++ )
 {
 double b[2];
 double a[2][2]; // 2*2
 
 int jPlus1 = ( j + 1 ) % nsRed;
 for( int k = 0; k < 2; k++ )
 {
 b[k] = red[2 * j + k] - blue[2 * i + k];
 // row k of a: a(k, 0), a(k, 1)
 a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
 a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
 }
 double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
 if( fabs( delta ) > 1.e-14 )
 {
 // not parallel
 double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
 double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
 if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
 {
 // the intersection is good
 for( int k = 0; k < 2; k++ )
 {
 points[2 * nPoints + k] =
 blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
 }
 markb[i] = 1; // so neighbor number i of blue will be considered too.
 markr[j] = 1; // this will be used in advancing red around blue quad
 nPoints++;
 }
 } // if the edges are too "parallel", skip them
 }
 }
 else // edge type is 1, so use 3d intersection
 {
 CartVect& C = bluec[i];
 CartVect& D = bluec[iPlus1];
 for( int j = 0; j < nsRed; j++ )
 {
 int jPlus1 = ( j + 1 ) % nsRed; // nsRed is just 4, forget about it, usually
 CartVect& A = redc[j];
 CartVect& B = redc[jPlus1];
 int np = 0;
 double E[9];
 intersect_great_circle_arc_with_clat_arc( A.array(), B.array(), C.array(), D.array(), R, E, np );
 if( np == 0 ) continue;
 if( np >= 2 )
 {
 std::cout << "intersection with 2 points :" << A << B << C << D << "\n";
 }
 for( int k = 0; k < np; k++ )
 {
 gnomonic_projection( CartVect( E + k * 3 ), R, plane, points[2 * nPoints],
 points[2 * nPoints + 1] );
 nPoints++;
 }
 markb[i] = 1; // so neighbor number i of blue will be considered too.
 markr[j] = 1; // this will be used in advancing red around blue quad
 }
 }
 }
 return MB_SUCCESS;
}

References moab::CartVect::array(), moab::E, gnomonic_projection(), intersect_great_circle_arc_with_clat_arc(), MB_SUCCESS, and moab::R.

Referenced by moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

enforce_convexity()

ErrorCode moab::IntxUtils::enforce_convexity ( Interface *  mb, EntityHandle  lset, int  my_rank = 0  )

static

Enforces convexity for a given set of polygons.

This function checks each polygon in the input set and computes the angles of each vertex. If a reflex angle is found, the polygon is broken into triangles and added back to the set. This process continues until all polygons in the set are convex.

Parameters

mb	The interface to the MOAB instance.
lset	The handle of the input set containing the polygons.
my_rank	The rank of the local process.

Returns

The error code indicating the success or failure of the operation.

Definition at line 1285 of file IntxUtils.cpp.

{
 Range inputRange;
 ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );MB_CHK_ERR( rval );
 
 Tag corrTag = 0;
 EntityHandle dumH = 0;
 rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
 if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
 
 Tag gidTag;
 rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
 
 std::vector< double > coords;
 coords.resize( 3 * MAXEDGES ); // at most 10 vertices per polygon
 // we should create a queue with new polygons that need processing for reflex angles
 // (obtuse)
 std::queue< EntityHandle > newPolys;
 int brokenPolys = 0;
 Range::iterator eit = inputRange.begin();
 while( eit != inputRange.end() || !newPolys.empty() )
 {
 EntityHandle eh;
 if( eit != inputRange.end() )
 {
 eh = *eit;
 ++eit;
 }
 else
 {
 eh = newPolys.front();
 newPolys.pop();
 }
 // get the nodes, then the coordinates
 const EntityHandle* verts;
 int num_nodes;
 rval = mb->get_connectivity( eh, verts, num_nodes );MB_CHK_ERR( rval );
 int nsides = num_nodes;
 // account for possible padded polygons
 while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
 nsides--;
 EntityHandle corrHandle = 0;
 if( corrTag )
 {
 rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );MB_CHK_ERR( rval );
 }
 int gid = 0;
 rval = mb->tag_get_data( gidTag, &eh, 1, &gid );MB_CHK_ERR( rval );
 coords.resize( 3 * nsides );
 if( nsides < 4 ) continue; // if already triangles, don't bother
 // get coordinates
 rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR( rval );
 // compute each angle
 bool alreadyBroken = false;
 
 for( int i = 0; i < nsides; i++ )
 {
 double* A = &coords[3 * i];
 double* B = &coords[3 * ( ( i + 1 ) % nsides )];
 double* C = &coords[3 * ( ( i + 2 ) % nsides )];
 double angle = IntxUtils::oriented_spherical_angle( A, B, C );
 if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
 {
 if( alreadyBroken )
 {
 mb->list_entities( &eh, 1 );
 mb->list_entities( verts, nsides );
 double* D = &coords[3 * ( ( i + 3 ) % nsides )];
 std::cout << "ABC: " << angle << " \n";
 std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
 std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
 std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
 std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
 
 return MB_FAILURE;
 }
 // the bad angle is at i+1;
 // create 1 triangle and one polygon; add the polygon to the input range, so
 // it will be processed too
 // also, add both to the set :) and remove the original polygon from the set
 // break the next triangle, even though not optimal
 // so create the triangle i+1, i+2, i+3; remove i+2 from original list
 // even though not optimal in general, it is good enough.
 EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
 verts[( i + 3 ) % nsides] };
 // create a polygon with num_nodes-1 vertices, and connectivity
 // verts[i+1], verts[i+3], (all except i+2)
 std::vector< EntityHandle > conn( nsides - 1 );
 for( int j = 1; j < nsides; j++ )
 {
 conn[j - 1] = verts[( i + j + 2 ) % nsides];
 }
 EntityHandle newElement;
 rval = mb->create_element( MBTRI, conn3, 3, newElement );MB_CHK_ERR( rval );
 
 rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
 if( corrTag )
 {
 rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
 }
 rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
 if( nsides == 4 )
 {
 // create another triangle
 rval = mb->create_element( MBTRI, &conn[0], 3, newElement );MB_CHK_ERR( rval );
 }
 else
 {
 // create another polygon, and add it to the inputRange
 rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );MB_CHK_ERR( rval );
 newPolys.push( newElement ); // because it has less number of edges, the
 // reverse should work to find it.
 }
 rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
 if( corrTag )
 {
 rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
 }
 rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
 rval = mb->remove_entities( lset, &eh, 1 );MB_CHK_ERR( rval );
 brokenPolys++;
 alreadyBroken = true; // get out of the loop, element is broken
 }
 }
 }
 if( brokenPolys > 0 )
 {
 std::cout << "on local process " << my_rank << ", " << brokenPolys
 << " concave polygons were decomposed in convex ones \n";
#ifdef VERBOSE
 std::stringstream fff;
 fff << "file_set" << mb->id_from_handle( lset ) << "rk_" << my_rank << ".h5m";
 rval = mb->write_file( fff.str().c_str(), 0, 0, &lset, 1 );MB_CHK_ERR( rval );
 std::cout << "wrote new file set: " << fff.str() << "\n";
#endif
 }
 return MB_SUCCESS;
}

References moab::angle(), moab::Range::begin(), CORRTAGNAME, moab::Range::end(), ErrorCode, MAXEDGES, mb, MB_CHK_ERR, MB_SUCCESS, MB_TAG_DENSE, MB_TAG_NOT_FOUND, MB_TYPE_HANDLE, MB_TYPE_INTEGER, MBPOLYGON, MBTRI, and oriented_spherical_angle().

Referenced by main().

fix_degenerate_quads()

ErrorCode moab::IntxUtils::fix_degenerate_quads ( Interface *  mb, EntityHandle  set  )

static

Definition at line 1426 of file IntxUtils.cpp.

{
 Range quads;
 ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
 Tag gid;
 gid = mb->globalId_tag();
 for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
 {
 EntityHandle quad = *qit;
 const EntityHandle* conn4 = NULL;
 int num_nodes = 0;
 rval = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
 for( int i = 0; i < num_nodes; i++ )
 {
 int next_node_index = ( i + 1 ) % num_nodes;
 if( conn4[i] == conn4[next_node_index] )
 {
 // form a triangle and delete the quad
 // first get the global id, to set it on triangle later
 int global_id = 0;
 rval = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
 int i2 = ( i + 2 ) % num_nodes;
 int i3 = ( i + 3 ) % num_nodes;
 EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
 EntityHandle tri;
 rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
 mb->add_entities( set, &tri, 1 );
 mb->remove_entities( set, &quad, 1 );
 mb->delete_entities( &quad, 1 );
 rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
 }
 }
 }
 return MB_SUCCESS;
}

References moab::Range::begin(), moab::Range::end(), ErrorCode, mb, MB_CHK_ERR, MB_SUCCESS, MBQUAD, and MBTRI.

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

global_gnomonic_projection()

ErrorCode moab::IntxUtils::global_gnomonic_projection ( Interface *  mb, EntityHandle  inSet, double  R, bool  centers_only, EntityHandle &  outSet  )

static

Definition at line 611 of file IntxUtils.cpp.

{
 std::string parTagName( "PARALLEL_PARTITION" );
 Tag part_tag;
 Range partSets;
 ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
 if( MB_SUCCESS == rval && part_tag != 0 )
 {
 rval = mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, NULL, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
 }
 rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
 // Get all entities of dimension 2
 Range inputRange; // get
 rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
 rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
 
 std::map< EntityHandle, int > partsAssign;
 std::map< int, EntityHandle > newPartSets;
 if( !partSets.empty() )
 {
 // get all cells, and assign parts
 for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
 {
 EntityHandle pSet = *setIt;
 Range ents;
 rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
 int val;
 rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
 // create a new set with the same part id tag, in the outSet
 EntityHandle newPartSet;
 rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
 rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
 newPartSets[val] = newPartSet;
 rval = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
 for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
 {
 partsAssign[*it] = val;
 }
 }
 }
 
 if( centers_only )
 {
 for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
 {
 CartVect center;
 EntityHandle cell = *it;
 rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
 int plane;
 decide_gnomonic_plane( center, plane );
 double c[3];
 c[2] = 0.;
 gnomonic_projection( center, R, plane, c[0], c[1] );
 
 gnomonic_unroll( c[0], c[1], R, plane );
 
 EntityHandle vertex;
 rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
 rval = mb->add_entities( outSet, &vertex, 1 );MB_CHK_ERR( rval );
 }
 }
 else
 {
 // distribute the cells to 6 planes, based on the center
 Range subranges[6];
 for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
 {
 CartVect center;
 EntityHandle cell = *it;
 rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
 int plane;
 decide_gnomonic_plane( center, plane );
 subranges[plane - 1].insert( cell );
 }
 for( int i = 1; i <= 6; i++ )
 {
 Range verts;
 rval = mb->get_connectivity( subranges[i - 1], verts );MB_CHK_ERR( rval );
 std::map< EntityHandle, EntityHandle > corr;
 for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
 {
 CartVect vect;
 EntityHandle v = *vt;
 rval = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
 double c[3];
 c[2] = 0.;
 gnomonic_projection( vect, R, i, c[0], c[1] );
 gnomonic_unroll( c[0], c[1], R, i );
 EntityHandle vertex;
 rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
 corr[v] = vertex; // for new connectivity
 }
 EntityHandle new_conn[20]; // max edges in 2d ?
 for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
 {
 EntityHandle eh = *eit;
 const EntityHandle* conn = NULL;
 int num_nodes;
 rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
 // build a new vertex array
 for( int j = 0; j < num_nodes; j++ )
 new_conn[j] = corr[conn[j]];
 EntityType type = mb->type_from_handle( eh );
 EntityHandle newCell;
 rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
 rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
 std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
 if( mit != partsAssign.end() )
 {
 int val = mit->second;
 rval = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
 }
 }
 }
 }
 
 return MB_SUCCESS;
}

References moab::CartVect::array(), moab::Range::begin(), center(), decide_gnomonic_plane(), moab::Range::empty(), moab::Range::end(), ErrorCode, gnomonic_projection(), gnomonic_unroll(), moab::Range::insert(), mb, MB_CHK_ERR, MB_SUCCESS, MBENTITYSET, MESHSET_SET, moab::R, ScaleToRadius(), and moab::Interface::UNION.

gnomonic_projection()

ErrorCode moab::IntxUtils::gnomonic_projection ( const CartVect &  pos, double  R, int  plane, double &  c1, double &  c2  )

static

Definition at line 456 of file IntxUtils.cpp.

{
 double alfa = 1.; // the new point will be on line alfa*pos
 
 switch( plane )
 {
 case 1: {
 // the plane with x = R; x>y, x>z
 // c1->y, c2->z
 alfa = R / pos[0];
 c1 = alfa * pos[1];
 c2 = alfa * pos[2];
 break;
 }
 case 2: {
 // y = R -> zx
 alfa = R / pos[1];
 c1 = alfa * pos[2];
 c2 = alfa * pos[0];
 break;
 }
 case 3: {
 // x=-R, -> yz
 alfa = -R / pos[0];
 c1 = -alfa * pos[1]; // the sign is to preserve orientation
 c2 = alfa * pos[2];
 break;
 }
 case 4: {
 // y = -R
 alfa = -R / pos[1];
 c1 = -alfa * pos[2]; // the sign is to preserve orientation
 c2 = alfa * pos[0];
 break;
 }
 case 5: {
 // z = -R
 alfa = -R / pos[2];
 c1 = -alfa * pos[0]; // the sign is to preserve orientation
 c2 = alfa * pos[1];
 break;
 }
 case 6: {
 alfa = R / pos[2];
 c1 = alfa * pos[0];
 c2 = alfa * pos[1];
 break;
 }
 default:
 return MB_FAILURE; // error
 }
 
 return MB_SUCCESS; // no error
}

References MB_SUCCESS, and moab::R.

Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), EdgeIntxRllCs(), global_gnomonic_projection(), moab::Intx2MeshOnSphere::setup_tgt_cell(), moab::IntxRllCssphere::setup_tgt_cell(), transform_coordinates(), and moab::BoundBox::update_box_spherical_elem().

gnomonic_unroll()

void moab::IntxUtils::gnomonic_unroll ( double &  c1, double &  c2, double  R, int  plane  )

static

Definition at line 574 of file IntxUtils.cpp.

{
 double tmp;
 switch( plane )
 {
 case 1:
 break; // nothing
 case 2: // rotate + 90
 tmp = c1;
 c1 = -c2;
 c2 = tmp;
 c1 = c1 + 2 * R;
 break;
 case 3:
 c1 = c1 + 4 * R;
 break;
 case 4: // rotate with -90 x-> -y; y -> x
 
 tmp = c1;
 c1 = c2;
 c2 = -tmp;
 c1 = c1 - 2 * R;
 break;
 case 5: // South Pole
 // rotate 180 then move to (-2, -2)
 c1 = -c1 - 2. * R;
 c2 = -c2 - 2. * R;
 break;
 ;
 case 6: // North Pole
 c1 = c1 - 2. * R;
 c2 = c2 + 2. * R;
 break;
 }
 return;
}

References moab::R.

Referenced by global_gnomonic_projection(), and transform_coordinates().

intersect_great_circle_arc_with_clat_arc()

ErrorCode moab::IntxUtils::intersect_great_circle_arc_with_clat_arc ( double *  A, double *  B, double *  C, double *  D, double  R, double *  E, int &  np  )

static

Definition at line 1605 of file IntxUtils.cpp.

{
 const double distTol = R * 1.e-6;
 const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
 np = 0; // number of points in intersection
 CartVect a( A ), b( B ), c( C ), d( D );
 // check input first
 double R2 = R * R;
 if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
 fabs( d.length_squared() - R2 ) >
 10 * Tolerance )
 return MB_FAILURE;
 
 if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
 if( ( c - d ).length_squared() < Tolerance ) // edges are too short
 return MB_FAILURE;
 
 // CD is the const latitude arc
 if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
 return MB_FAILURE;
 
 if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
 
 // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
 // great circle arc AB normal to the AB circle:
 CartVect n1 = a * b; // the normal to the great circle arc (circle)
 // solve the system of equations:
 /*
 * n1%(x, y, z) = 0 // on the great circle
 * z = C[2];
 * x^2+y^2+z^2 = R^2
 */
 double z = C[2];
 if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
 {
 // it is the Equator; check if the const lat edge is Equator too
 if( fabs( C[2] ) > distTol )
 {
 return MB_FAILURE; // no intx, too far from Eq
 }
 else
 {
 // all points are on the equator
 //
 CartVect cd = c * d;
 // by convention, c<d, positive is from c to d
 // is a or b between c , d?
 CartVect ca = c * a;
 CartVect ad = a * d;
 CartVect cb = c * b;
 CartVect bd = b * d;
 bool agtc = ( ca % cd >= -Tolerance ); // a>c?
 bool dgta = ( ad % cd >= -Tolerance ); // d>a?
 bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
 bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
 if( agtc )
 {
 if( dgta )
 {
 // a is for sure a point
 E[0] = a[0];
 E[1] = a[1];
 E[2] = a[2];
 np++;
 if( bgtc )
 {
 if( dgtb )
 {
 // b is also in between c and d
 E[3] = b[0];
 E[4] = b[1];
 E[5] = b[2];
 np++;
 }
 else
 {
 // then order is c a d b, intx is ad
 E[3] = d[0];
 E[4] = d[1];
 E[5] = d[2];
 np++;
 }
 }
 else
 {
 // b is less than c, so b c a d, intx is ac
 E[3] = c[0];
 E[4] = c[1];
 E[5] = c[2];
 np++; // what if E[0] is E[3]?
 }
 }
 else // c < d < a
 {
 if( dgtb ) // d is for sure in
 {
 E[0] = d[0];
 E[1] = d[1];
 E[2] = d[2];
 np++;
 if( bgtc ) // c<b<d<a
 {
 // another point is b
 E[3] = b[0];
 E[4] = b[1];
 E[5] = b[2];
 np++;
 }
 else // b<c<d<a
 {
 // another point is c
 E[3] = c[0];
 E[4] = c[1];
 E[5] = c[2];
 np++;
 }
 }
 else
 {
 // nothing, order is c, d < a, b
 }
 }
 }
 else // a < c < d
 {
 if( bgtc )
 {
 // c is for sure in
 E[0] = c[0];
 E[1] = c[1];
 E[2] = c[2];
 np++;
 if( dgtb )
 {
 // a < c < b < d; second point is b
 E[3] = b[0];
 E[4] = b[1];
 E[5] = b[2];
 np++;
 }
 else
 {
 // a < c < d < b; second point is d
 E[3] = d[0];
 E[4] = d[1];
 E[5] = d[2];
 np++;
 }
 }
 else // a, b < c < d
 {
 // nothing
 }
 }
 }
 // for the 2 points selected, see if it is only one?
 // no problem, maybe it will be collapsed later anyway
 if( np > 0 ) return MB_SUCCESS;
 return MB_FAILURE; // no intersection
 }
 {
 if( fabs( n1[0] ) <= fabs( n1[1] ) )
 {
 // resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
 // (u+v*x)^2+x^2=R2-z^2
 // (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
 // delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
 // x1,2 =
 double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
 double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
 double delta = b1 * b1 - 4 * a1 * c1;
 if( delta < -Tolerance ) return MB_FAILURE; // no intersection
 if( delta > Tolerance ) // 2 solutions possible
 {
 double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
 double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
 double y1 = u + v * x1;
 double y2 = u + v * x2;
 if( verify( a, b, c, d, x1, y1, z ) )
 {
 E[0] = x1;
 E[1] = y1;
 E[2] = z;
 np++;
 }
 if( verify( a, b, c, d, x2, y2, z ) )
 {
 E[3 * np + 0] = x2;
 E[3 * np + 1] = y2;
 E[3 * np + 2] = z;
 np++;
 }
 }
 else
 {
 // one solution
 double x1 = -b1 / 2 / a1;
 double y1 = u + v * x1;
 if( verify( a, b, c, d, x1, y1, z ) )
 {
 E[0] = x1;
 E[1] = y1;
 E[2] = z;
 np++;
 }
 }
 }
 else
 {
 // resolve eq in y, reverse
 // n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
 // (u+v*y)^2+y^2 -R2+z^2 =0
 // (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
 //
 // x1,2 =
 double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
 double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
 double delta = b1 * b1 - 4 * a1 * c1;
 if( delta < -Tolerance ) return MB_FAILURE; // no intersection
 if( delta > Tolerance ) // 2 solutions possible
 {
 double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
 double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
 double x1 = u + v * y1;
 double x2 = u + v * y2;
 if( verify( a, b, c, d, x1, y1, z ) )
 {
 E[0] = x1;
 E[1] = y1;
 E[2] = z;
 np++;
 }
 if( verify( a, b, c, d, x2, y2, z ) )
 {
 E[3 * np + 0] = x2;
 E[3 * np + 1] = y2;
 E[3 * np + 2] = z;
 np++;
 }
 }
 else
 {
 // one solution
 double y1 = -b1 / 2 / a1;
 double x1 = u + v * y1;
 if( verify( a, b, c, d, x1, y1, z ) )
 {
 E[0] = x1;
 E[1] = y1;
 E[2] = z;
 np++;
 }
 }
 }
 }
 
 if( np <= 0 ) return MB_FAILURE;
 return MB_SUCCESS;
}

References moab::E, moab::CartVect::length_squared(), length_squared(), MB_SUCCESS, moab::R, and moab::verify().

Referenced by EdgeIntxRllCs().

intersect_great_circle_arcs()

ErrorCode moab::IntxUtils::intersect_great_circle_arcs ( double *  A, double *  B, double *  C, double *  D, double  R, double *  E  )

static

Definition at line 1518 of file IntxUtils.cpp.

{
 // first verify A, B, C, D are on the same sphere
 double R2 = R * R;
 const double Tolerance = 1.e-12 * R2;
 
 CartVect a( A ), b( B ), c( C ), d( D );
 
 if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
 fabs( d.length_squared() - R2 ) >
 10 * Tolerance )
 return MB_FAILURE;
 
 CartVect n1 = a * b;
 if( n1.length_squared() < Tolerance ) return MB_FAILURE;
 
 CartVect n2 = c * d;
 if( n2.length_squared() < Tolerance ) return MB_FAILURE;
 CartVect n3 = n1 * n2;
 n3.normalize();
 
 n3 = R * n3;
 // the intersection is either n3 or -n3
 CartVect n4 = a * n3, n5 = n3 * b;
 if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
 {
 // n3 is good for ab, see if it is good for cd
 n4 = c * n3;
 n5 = n3 * d;
 if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
 {
 E[0] = n3[0];
 E[1] = n3[1];
 E[2] = n3[2];
 }
 else
 return MB_FAILURE;
 }
 else
 {
 // try -n3
 n3 = -n3;
 n4 = a * n3, n5 = n3 * b;
 if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
 {
 // n3 is good for ab, see if it is good for cd
 n4 = c * n3;
 n5 = n3 * d;
 if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
 {
 E[0] = n3[0];
 E[1] = n3[1];
 E[2] = n3[2];
 }
 else
 return MB_FAILURE;
 }
 else
 return MB_FAILURE;
 }
 
 return MB_SUCCESS;
}

References moab::E, moab::CartVect::length_squared(), MB_SUCCESS, moab::CartVect::normalize(), and moab::R.

oriented_spherical_angle()

double moab::IntxUtils::oriented_spherical_angle ( const double *  A, const double *  B, const double *  C  )

static

Definition at line 876 of file IntxUtils.cpp.

{
 // assume the same radius, sphere at origin
 CartVect a( A ), b( B ), c( C );
 CartVect normalOAB = a * b;
 CartVect normalOCB = c * b;
 CartVect orient = ( c - b ) * ( a - b );
 double ang = angle( normalOAB, normalOCB ); // this is between 0 and M_PI
 if( ang != ang )
 {
 // signal of a nan
 std::cout << a << " " << b << " " << c << "\n";
 std::cout << ang << "\n";
 }
 if( orient % b < 0 ) return ( 2 * M_PI - ang ); // the other angle, supplement
 
 return ang;
}

References moab::angle().

Referenced by moab::IntxAreaUtils::area_spherical_polygon_girard(), and enforce_convexity().

remove_duplicate_vertices()

ErrorCode moab::IntxUtils::remove_duplicate_vertices ( Interface *  mb, EntityHandle  file_set, double  merge_tol, std::vector< Tag > &  tagList  )

static

Definition at line 2062 of file IntxUtils.cpp.

{
 Range verts;
 ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
 rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
 
 MergeMesh mm( mb );
 
 // remove the vertices from the set, before merging
 
 rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
 
 // now correct vertices that are repeated in polygons
 rval = remove_padded_vertices( mb, file_set, tagList );
 return MB_SUCCESS;
}

References ErrorCode, mb, MB_CHK_ERR, MB_SUCCESS, moab::MergeMesh::merge_all(), and remove_padded_vertices().

remove_padded_vertices()

ErrorCode moab::IntxUtils::remove_padded_vertices ( Interface *  mb, EntityHandle  file_set, std::vector< Tag > &  tagList  )

static

Definition at line 2082 of file IntxUtils.cpp.

{
 
 // now correct vertices that are repeated in polygons
 Range cells;
 ErrorCode rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
 
 Range verts;
 rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
 
 Range modifiedCells; // will be deleted at the end; keep the gid
 Range newCells;
 
 for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
 {
 EntityHandle cell = *cit;
 const EntityHandle* connec = NULL;
 int num_verts = 0;
 rval = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
 
 std::vector< EntityHandle > newConnec;
 newConnec.push_back( connec[0] ); // at least one vertex
 int index = 0;
 int new_size = 1;
 while( index < num_verts - 2 )
 {
 int next_index = ( index + 1 );
 if( connec[next_index] != newConnec[new_size - 1] )
 {
 newConnec.push_back( connec[next_index] );
 new_size++;
 }
 index++;
 }
 // add the last one only if different from previous and first node
 if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
 {
 newConnec.push_back( connec[num_verts - 1] );
 new_size++;
 }
 if( new_size < num_verts )
 {
 // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
 modifiedCells.insert( cell );
 // create a new cell with type triangle, quad or polygon
 EntityType type = MBTRI;
 if( new_size == 3 )
 type = MBTRI;
 else if( new_size == 4 )
 type = MBQUAD;
 else if( new_size > 4 )
 type = MBPOLYGON;
 
 // create new cell
 EntityHandle newCell;
 rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
 // set the old id to the new element
 newCells.insert( newCell );
 double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
 for( size_t i = 0; i < tagList.size(); i++ )
 {
 rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
 rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
 }
 }
 }
 
 rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
 rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
 rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
 rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
 
 return MB_SUCCESS;
}

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Range::insert(), mb, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MBPOLYGON, MBQUAD, and MBTRI.

Referenced by moab::NCHelperDomain::create_mesh(), moab::NCHelperScrip::create_mesh(), and remove_duplicate_vertices().

reverse_gnomonic_projection()

ErrorCode moab::IntxUtils::reverse_gnomonic_projection ( const double &  c1, const double &  c2, double  R, int  plane, CartVect &  pos  )

static

Definition at line 512 of file IntxUtils.cpp.

{
 
 // the new point will be on line beta*pos
 double len = sqrt( c1 * c1 + c2 * c2 + R * R );
 double beta = R / len; // it is less than 1, in general
 
 switch( plane )
 {
 case 1: {
 // the plane with x = R; x>y, x>z
 // c1->y, c2->z
 pos[0] = beta * R;
 pos[1] = c1 * beta;
 pos[2] = c2 * beta;
 break;
 }
 case 2: {
 // y = R -> zx
 pos[1] = R * beta;
 pos[2] = c1 * beta;
 pos[0] = c2 * beta;
 break;
 }
 case 3: {
 // x=-R, -> yz
 pos[0] = -R * beta;
 pos[1] = -c1 * beta; // the sign is to preserve orientation
 pos[2] = c2 * beta;
 break;
 }
 case 4: {
 // y = -R
 pos[1] = -R * beta;
 pos[2] = -c1 * beta; // the sign is to preserve orientation
 pos[0] = c2 * beta;
 break;
 }
 case 5: {
 // z = -R
 pos[2] = -R * beta;
 pos[0] = -c1 * beta; // the sign is to preserve orientation
 pos[1] = c2 * beta;
 break;
 }
 case 6: {
 pos[2] = R * beta;
 pos[0] = c1 * beta;
 pos[1] = c2 * beta;
 break;
 }
 default:
 return MB_FAILURE; // error
 }
 
 return MB_SUCCESS; // no error
}

References MB_SUCCESS, and moab::R.

Referenced by moab::Intx2MeshOnSphere::findNodes(), moab::IntxRllCssphere::findNodes(), and moab::BoundBox::update_box_spherical_elem().

ScaleToRadius()

ErrorCode moab::IntxUtils::ScaleToRadius ( Interface *  mb, EntityHandle  set, double  R  )

static

Definition at line 833 of file IntxUtils.cpp.

{
 Range nodes;
 ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true ); // recursive
 if( rval != moab::MB_SUCCESS ) return rval;
 
 // one by one, get the node and project it on the sphere, with a radius given
 // the center of the sphere is at 0,0,0
 for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
 {
 EntityHandle nd = *nit;
 CartVect pos;
 rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
 if( rval != moab::MB_SUCCESS ) return rval;
 double len = pos.length();
 if( len == 0. ) return MB_FAILURE;
 pos = R / len * pos;
 rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
 if( rval != moab::MB_SUCCESS ) return rval;
 }
 return MB_SUCCESS;
}

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::CartVect::length(), mb, MB_SUCCESS, MBVERTEX, and moab::R.

Referenced by CreateTempestMesh(), global_gnomonic_projection(), and main().

SortAndRemoveDoubles2()

int moab::IntxUtils::SortAndRemoveDoubles2 ( double *  P, int &  nP, double  epsilon_1  )

static

Sorts and removes duplicate points in the given array.

Note: nP might be modified too, we will remove duplicates if found

Parameters

P	The array of points to be sorted and checked for duplicates.
nP	The number of points in P.
epsilon_1	The epsilon value for distance comparison.

Returns

0 if successful.

Definition at line 130 of file IntxUtils.cpp.

{
 if( nP < 2 ) return 0; // nothing to do
 
 // center of gravity for the points
 double c[2] = { 0., 0. };
 int k = 0;
 for( k = 0; k < nP; k++ )
 {
 c[0] += P[2 * k];
 c[1] += P[2 * k + 1];
 }
 c[0] /= nP;
 c[1] /= nP;
 
 // how many? we dimensioned P at MAXEDGES*10; so we imply we could have at most 5*MAXEDGES
 // intersection points
 struct angleAndIndex pairAngleIndex[5 * MAXEDGES];
 
 for( k = 0; k < nP; k++ )
 {
 double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
 if( x != 0. || y != 0. )
 {
 pairAngleIndex[k].angle = atan2( y, x );
 }
 else
 {
 pairAngleIndex[k].angle = 0;
 // it would mean that the cells are touching at a vertex
 }
 pairAngleIndex[k].index = k;
 }
 
 // this should be faster than the bubble sort we had before
 std::sort( pairAngleIndex, pairAngleIndex + nP, angleCompare );
 // copy now to a new double array
 double PCopy[10 * MAXEDGES]; // the same dimension as P; very conservative, but faster than
 // reallocate for a vector
 for( k = 0; k < nP; k++ ) // index will show where it should go now;
 {
 int ck = pairAngleIndex[k].index;
 PCopy[2 * k] = P[2 * ck];
 PCopy[2 * k + 1] = P[2 * ck + 1];
 }
 // now copy from PCopy over original P location
 std::copy( PCopy, PCopy + 2 * nP, P );
 
 // eliminate duplicates, finally
 
 int i = 0, j = 1; // the next one; j may advance faster than i
 // check the unit
 // double epsilon_1 = 1.e-5; // these are cm; 2 points are the same if the distance is less
 // than 1.e-5 cm
 while( j < nP )
 {
 double d2 = dist2( &P[2 * i], &P[2 * j] );
 if( d2 > epsilon_1 )
 {
 i++;
 P[2 * i] = P[2 * j];
 P[2 * i + 1] = P[2 * j + 1];
 }
 j++;
 }
 // test also the last point with the first one (index 0)
 // the first one could be at -PI; last one could be at +PI, according to atan2 span
 
 double d2 = dist2( P, &P[2 * i] ); // check the first and last points (ordered from -pi to +pi)
 if( d2 > epsilon_1 )
 {
 nP = i + 1;
 }
 else
 nP = i; // effectively delete the last point (that would have been the same with first)
 if( nP == 0 ) nP = 1; // we should be left with at least one point we already tested if nP is 0 originally
 return 0;
}

References moab::angleAndIndex::angle, moab::angleCompare(), dist2(), moab::angleAndIndex::index, and MAXEDGES.

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc(), and moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc().

spherical_to_cart()

CartVect moab::IntxUtils::spherical_to_cart ( IntxUtils::SphereCoords &  sc )

static

Definition at line 824 of file IntxUtils.cpp.

{
 CartVect res;
 res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
 res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
 res[2] = sc.R * sin( sc.lat ); // z
 return res;
}

References moab::IntxUtils::SphereCoords::lat, moab::IntxUtils::SphereCoords::lon, and moab::IntxUtils::SphereCoords::R.

Referenced by main().

transform_coordinates()

void moab::IntxUtils::transform_coordinates ( double *  avg_position, int  projection_type  )

static

Definition at line 733 of file IntxUtils.cpp.

{
 if( projection_type == 1 )
 {
 double R =
 avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
 R = sqrt( R );
 double lat = asin( avg_position[2] / R );
 double lon = atan2( avg_position[1], avg_position[0] );
 avg_position[0] = lon;
 avg_position[1] = lat;
 avg_position[2] = R;
 }
 else if( projection_type == 2 ) // gnomonic projection
 {
 CartVect pos( avg_position );
 int gplane;
 IntxUtils::decide_gnomonic_plane( pos, gplane );
 
 IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
 avg_position[2] = 0;
 IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
 }
}

References decide_gnomonic_plane(), gnomonic_projection(), gnomonic_unroll(), and moab::R.

---

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

• IntxUtils.hpp
• IntxUtils.cpp